tutorial/tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>SpaGCN Tutorial</center></h1>\n",
    "\n",
    "\n",
    "<center>Author: Jian Hu,*, Xiangjie Li, Kyle Coleman, Amelia Schroeder, Nan Ma, David J. Irwin, Edward B. Lee, Russell T. Shinohara, Mingyao Li*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outline\n",
    "1. Installation\n",
    "2. Import modules\n",
    "3. Read in data\n",
    "4. Integrate gene expression and histology into a Graph\n",
    "5. Spatial domain detection using SpaGCN\n",
    "6. Identify SVGs\n",
    "7. Identify Meta Gene\n",
    "8. Multiple tissue sections analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Installation\n",
    "To install SpaGCN package you must make sure that your python version is over 3.5.=. If you don’t know the version of python you can check it by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'3.8.8'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import platform\n",
    "platform.python_version()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Because SpaGCN pends on pytorch, you should make sure torch is correctly installed.\n",
    "<br>\n",
    "Now you can install the current release of SpaGCN by the following three ways:\n",
    "#### 1.1 PyPI: Directly install the package from PyPI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip3 install SpaGCN\n",
    "#Note: you need to make sure that the pip is for python3，or we should install SpaGCN by\n",
    "python3 -m pip install SpaGCN\n",
    "pip3 install SpaGCN\n",
    "#If you do not have permission (when you get a permission denied error), you should install SpaGCN by\n",
    "pip3 install --user SpaGCN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2 Github\n",
    "Download the package from Github and install it locally:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "git clone https://github.com/jianhuupenn/SpaGCN\n",
    "cd SpaGCN/SpaGCN_package/\n",
    "python3 setup.py install --user"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.3 Anaconda\n",
    "If you do not have Python3.5 or Python3.6 installed, consider installing Anaconda (see Installing Anaconda). After installing Anaconda, you can create a new environment, for example, SpaGCN (you can change to any name you like)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create an environment called SpaGCN\n",
    "conda create -n SpaGCN python=3.7.9\n",
    "#activate your environment \n",
    "conda activate SpaGCN\n",
    "git clone https://github.com/jianhuupenn/SpaGCN\n",
    "cd SpaGCN/SpaGCN_package/\n",
    "python3 setup.py build\n",
    "python3 setup.py install\n",
    "conda deactivate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Import python modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os,csv,re\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scanpy as sc\n",
    "import math\n",
    "import SpaGCN as spg\n",
    "from scipy.sparse import issparse\n",
    "import random, torch\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import matplotlib.colors as clr\n",
    "import matplotlib.pyplot as plt\n",
    "import SpaGCN as spg\n",
    "#In order to read in image data, we need to install some package. Here we recommend package \"opencv\"\n",
    "#inatll opencv in python\n",
    "#!pip3 install opencv-python\n",
    "import cv2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.2.0'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spg.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Read in data\n",
    "The current version of SpaGCN requres three input data: \n",
    "<br>\n",
    "1. The gene expression matrix(n by k): expression_matrix.h5;\n",
    "<br>\n",
    "2. Spatial coordinateds of samplespositions.txt;\n",
    "<br>\n",
    "3. Histology image(optional): histology.tif, can be tif or png or jepg.\n",
    "<br>\n",
    "The gene expreesion data can be stored as an AnnData object. AnnData stores a data matrix .X together with annotations of observations .obs, variables .var and unstructured annotations .uns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "#Read original 10x_h5 data and save it to h5ad\n",
    "from scanpy import read_10x_h5\n",
    "adata = read_10x_h5(\"../tutorial/data/151673/expression_matrix.h5\")\n",
    "spatial=pd.read_csv(\"../tutorial/data/151673/positions.txt\",sep=\",\",header=None,na_filter=False,index_col=0) \n",
    "adata.obs[\"x1\"]=spatial[1]\n",
    "adata.obs[\"x2\"]=spatial[2]\n",
    "adata.obs[\"x3\"]=spatial[3]\n",
    "adata.obs[\"x4\"]=spatial[4]\n",
    "adata.obs[\"x5\"]=spatial[5]\n",
    "adata.obs[\"x_array\"]=adata.obs[\"x2\"]\n",
    "adata.obs[\"y_array\"]=adata.obs[\"x3\"]\n",
    "adata.obs[\"x_pixel\"]=adata.obs[\"x4\"]\n",
    "adata.obs[\"y_pixel\"]=adata.obs[\"x5\"]\n",
    "#Select captured samples\n",
    "adata=adata[adata.obs[\"x1\"]==1]\n",
    "adata.var_names=[i.upper() for i in list(adata.var_names)]\n",
    "adata.var[\"genename\"]=adata.var.index.astype(\"str\")\n",
    "adata.write_h5ad(\"../tutorial/data/151673/sample_data.h5ad\")\n",
    "\"\"\"\n",
    "#Read in gene expression and spatial location\n",
    "adata=sc.read(\"../tutorial/data/151673/sample_data.h5ad\")\n",
    "#Read in hitology image\n",
    "img=cv2.imread(\"../tutorial/data/151673/histology.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Integrate gene expression and histology into a Graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Set coordinates\n",
    "x_array=adata.obs[\"x_array\"].tolist()\n",
    "y_array=adata.obs[\"y_array\"].tolist()\n",
    "x_pixel=adata.obs[\"x_pixel\"].tolist()\n",
    "y_pixel=adata.obs[\"y_pixel\"].tolist()\n",
    "\n",
    "#Test coordinates on the image\n",
    "img_new=img.copy()\n",
    "for i in range(len(x_pixel)):\n",
    "    x=x_pixel[i]\n",
    "    y=y_pixel[i]\n",
    "    img_new[int(x-20):int(x+20), int(y-20):int(y+20),:]=0\n",
    "\n",
    "cv2.imwrite('./sample_results/151673_map.jpg', img_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The ‘s’ parameter determines the weight given to histology when calculating Euclidean distance between every two spots. ‘s = 1’ means that the histology pixel intensity value has the same scale variance as the (x,y) coordinates, whereas higher value of ‘s’ indicates higher scale variance, hence, higher weight to histology, when calculating the Euclidean distance. \n",
    "\n",
    "- The \"b\"parameter determines the area of each spot when extracting color intensity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculateing adj matrix using histology image...\n",
      "Var of c0,c1,c2 =  33.30687202862215 174.55510595352243 46.84205750749746\n",
      "Var of x,y,z =  5606737.526317932 4468793.817921193 5606737.526317932\n"
     ]
    }
   ],
   "source": [
    "#Calculate adjacent matrix\n",
    "s=1\n",
    "b=49\n",
    "adj=spg.calculate_adj_matrix(x=x_pixel,y=y_pixel, x_pixel=x_pixel, y_pixel=y_pixel, image=img, beta=b, alpha=s, histology=True)\n",
    "#If histlogy image is not available, SpaGCN can calculate the adjacent matrix using the fnction below\n",
    "#adj=calculate_adj_matrix(x=x_pixel,y=y_pixel, histology=False)\n",
    "np.savetxt('./data/151673/adj.csv', adj, delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Spatial domain detection using SpaGCN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1 Expression data preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n"
     ]
    }
   ],
   "source": [
    "adata=sc.read(\"./data/151673/sample_data.h5ad\")\n",
    "adj=np.loadtxt('./data/151673/adj.csv', delimiter=',')\n",
    "adata.var_names_make_unique()\n",
    "spg.prefilter_genes(adata,min_cells=3) # avoiding all genes are zeros\n",
    "spg.prefilter_specialgenes(adata)\n",
    "#Normalize and take log for UMI\n",
    "sc.pp.normalize_per_cell(adata)\n",
    "sc.pp.log1p(adata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2 Set hyper-parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- p: Percentage of total expression contributed by neighborhoods.\n",
    "- l: Parameter to control p."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run 1: l [0.01, 1000], p [0.0, 153.8820492650696]\n",
      "Run 2: l [0.01, 500.005], p [0.0, 28.015447343094223]\n",
      "Run 3: l [0.01, 250.0075], p [0.0, 4.240330523308446]\n",
      "Run 4: l [0.01, 125.00874999999999], p [0.0, 0.5157276735032843]\n",
      "Run 5: l [62.509375, 125.00874999999999], p [0.028496868560644373, 0.5157276735032843]\n",
      "Run 6: l [93.7590625, 125.00874999999999], p [0.18753135107474428, 0.5157276735032843]\n",
      "Run 7: l [109.38390625, 125.00874999999999], p [0.32801349789332424, 0.5157276735032843]\n",
      "Run 8: l [117.196328125, 125.00874999999999], p [0.4156469508032292, 0.5157276735032843]\n",
      "Run 9: l [121.1025390625, 125.00874999999999], p [0.4640926787304587, 0.5157276735032843]\n",
      "Run 10: l [123.05564453125, 125.00874999999999], p [0.48950676051756026, 0.5157276735032843]\n",
      "recommended l =  124.032197265625\n"
     ]
    }
   ],
   "source": [
    "p=0.5 \n",
    "#Find the l value given p\n",
    "l=spg.search_l(p, adj, start=0.01, end=1000, tol=0.01, max_run=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- n_clusters: Number of spatial domains wanted.\n",
    "- res: Resolution in the initial Louvain's Clustering methods. If the number of clusters is known, we can use the spg.search_res() fnction to search for suitable resolution(optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start at res =  0.7 step =  0.1\n",
      "Initializing cluster centers with louvain, resolution =  0.7\n",
      "Epoch  0\n",
      "Epoch  10\n",
      "Res =  0.7 Num of clusters =  7\n",
      "recommended res =  0.7\n"
     ]
    }
   ],
   "source": [
    "#If the number of clusters known, we can use the spg.search_res() fnction to search for suitable resolution(optional)\n",
    "#For this toy data, we set the number of clusters=7 since this tissue has 7 layers\n",
    "n_clusters=7\n",
    "#Set seed\n",
    "r_seed=t_seed=n_seed=100\n",
    "#Seaech for suitable resolution\n",
    "res=spg.search_res(adata, adj, l, n_clusters, start=0.7, step=0.1, tol=5e-3, lr=0.05, max_epochs=20, r_seed=r_seed, t_seed=t_seed, n_seed=n_seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.3 Run SpaGCN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing cluster centers with louvain, resolution =  0.7\n",
      "Epoch  0\n",
      "Epoch  10\n",
      "Epoch  20\n",
      "Epoch  30\n",
      "Epoch  40\n",
      "Epoch  50\n",
      "Epoch  60\n",
      "Epoch  70\n",
      "Epoch  80\n",
      "Epoch  90\n",
      "Epoch  100\n",
      "Epoch  110\n",
      "Epoch  120\n",
      "Epoch  130\n",
      "Epoch  140\n",
      "Epoch  150\n",
      "Epoch  160\n",
      "Epoch  170\n",
      "Epoch  180\n",
      "Epoch  190\n",
      "Calculateing adj matrix using xy only...\n"
     ]
    }
   ],
   "source": [
    "clf=spg.SpaGCN()\n",
    "clf.set_l(l)\n",
    "#Set seed\n",
    "random.seed(r_seed)\n",
    "torch.manual_seed(t_seed)\n",
    "np.random.seed(n_seed)\n",
    "#Run\n",
    "clf.train(adata,adj,init_spa=True,init=\"louvain\",res=res, tol=5e-3, lr=0.05, max_epochs=200)\n",
    "y_pred, prob=clf.predict()\n",
    "adata.obs[\"pred\"]= y_pred\n",
    "adata.obs[\"pred\"]=adata.obs[\"pred\"].astype('category')\n",
    "#Do cluster refinement(optional)\n",
    "#shape=\"hexagon\" for Visium data, \"square\" for ST data.\n",
    "adj_2d=spg.calculate_adj_matrix(x=x_array,y=y_array, histology=False)\n",
    "refined_pred=spg.refine(sample_id=adata.obs.index.tolist(), pred=adata.obs[\"pred\"].tolist(), dis=adj_2d, shape=\"hexagon\")\n",
    "adata.obs[\"refined_pred\"]=refined_pred\n",
    "adata.obs[\"refined_pred\"]=adata.obs[\"refined_pred\"].astype('category')\n",
    "#Save results\n",
    "adata.write_h5ad(\"./sample_results/results.h5ad\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.4 Plot spatial domains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata=sc.read(\"./sample_results/results.h5ad\")\n",
    "#Set colors used\n",
    "plot_color=[\"#F56867\",\"#FEB915\",\"#C798EE\",\"#59BE86\",\"#7495D3\",\"#D1D1D1\",\"#6D1A9C\",\"#15821E\",\"#3A84E6\",\"#997273\",\"#787878\",\"#DB4C6C\",\"#9E7A7A\",\"#554236\",\"#AF5F3C\",\"#93796C\",\"#F9BD3F\",\"#DAB370\",\"#877F6C\",\"#268785\"]\n",
    "#Plot spatial domains\n",
    "domains=\"pred\"\n",
    "num_celltype=len(adata.obs[domains].unique())\n",
    "adata.uns[domains+\"_colors\"]=list(plot_color[:num_celltype])\n",
    "ax=sc.pl.scatter(adata,alpha=1,x=\"y_pixel\",y=\"x_pixel\",color=domains,title=domains,color_map=plot_color,show=False,size=100000/adata.shape[0])\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.axes.invert_yaxis()\n",
    "plt.savefig(\"./sample_results/pred.png\", dpi=600)\n",
    "plt.close()\n",
    "\n",
    "#Plot refined spatial domains\n",
    "domains=\"refined_pred\"\n",
    "num_celltype=len(adata.obs[domains].unique())\n",
    "adata.uns[domains+\"_colors\"]=list(plot_color[:num_celltype])\n",
    "ax=sc.pl.scatter(adata,alpha=1,x=\"y_pixel\",y=\"x_pixel\",color=domains,title=domains,color_map=plot_color,show=False,size=100000/adata.shape[0])\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.axes.invert_yaxis()\n",
    "plt.savefig(\"./sample_results/refined_pred.png\", dpi=600)\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Spatial Domains**![](./sample_results/pred.png) **Refined Spatial Domains**![](./sample_results/refined_pred.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Identify SVGs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n"
     ]
    }
   ],
   "source": [
    "#Read in raw data\n",
    "raw=sc.read(\"../tutorial/data/151673/sample_data.h5ad\")\n",
    "raw.var_names_make_unique()\n",
    "raw.obs[\"pred\"]=adata.obs[\"pred\"].astype('category')\n",
    "raw.obs[\"x_array\"]=raw.obs[\"x2\"]\n",
    "raw.obs[\"y_array\"]=raw.obs[\"x3\"]\n",
    "raw.obs[\"x_pixel\"]=raw.obs[\"x4\"]\n",
    "raw.obs[\"y_pixel\"]=raw.obs[\"x5\"]\n",
    "#Convert sparse matrix to non-sparse\n",
    "raw.X=(raw.X.A if issparse(raw.X) else raw.X)\n",
    "raw.raw=raw\n",
    "sc.pp.log1p(raw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- target: Target domain to identify SVGs.\n",
    "- min_in_group_fraction: Minium in-group expression fraction.\n",
    "- min_in_out_group_ratio: Miniumn (in-group expression fraction) / (out-group expression fraction).\n",
    "- min_fold_change: Miniumn (in-group expression) / (out-group expression).\n",
    "- r: Radius to detect a spot's neighboring spots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculateing adj matrix using xy only...\n",
      "Calculateing adj matrix using xy only...\n",
      "Calculateing adj matrix using xy only...\n",
      "Run 1: radius [1.4142135381698608, 16.970561981201172], num_nbr [1.0, 321.149863760218]\n",
      "Calculateing adj matrix using xy only...\n",
      "Run 2: radius [1.4142135381698608, 9.192387759685516], num_nbr [1.0, 117.8283378746594]\n",
      "Calculateing adj matrix using xy only...\n",
      "Run 3: radius [1.4142135381698608, 5.303300648927689], num_nbr [1.0, 41.85013623978202]\n",
      "Calculateing adj matrix using xy only...\n",
      "Run 4: radius [1.4142135381698608, 3.3587570935487747], num_nbr [1.0, 20.04632152588556]\n",
      "Calculateing adj matrix using xy only...\n",
      "Run 5: radius [2.386485315859318, 3.3587570935487747], num_nbr [8.7574931880109, 20.04632152588556]\n",
      "Calculateing adj matrix using xy only...\n",
      "recommended radius =  2.8726212047040462 num_nbr=12.524523160762943\n",
      "radius= 2.8726212047040462 average number of neighbors for each spot is 12.524523160762943\n",
      " Cluster 0 has neighbors:\n",
      "Dmain  3 :  863\n",
      "Dmain  2 :  517\n",
      "SVGs for domain  0 : ['CAMK2N1', 'ENC1', 'GPM6A', 'ARPP19', 'HPCAL1']\n"
     ]
    }
   ],
   "source": [
    "#Use domain 0 as an example\n",
    "target=0\n",
    "#Set filtering criterials\n",
    "min_in_group_fraction=0.8\n",
    "min_in_out_group_ratio=1\n",
    "min_fold_change=1.5\n",
    "#Search radius such that each spot in the target domain has approximately 10 neighbors on average\n",
    "adj_2d=spg.calculate_adj_matrix(x=x_array, y=y_array, histology=False)\n",
    "start, end= np.quantile(adj_2d[adj_2d!=0],q=0.001), np.quantile(adj_2d[adj_2d!=0],q=0.1)\n",
    "r=spg.search_radius(target_cluster=target, cell_id=adata.obs.index.tolist(), x=x_array, y=y_array, pred=adata.obs[\"pred\"].tolist(), start=start, end=end, num_min=10, num_max=14,  max_run=100)\n",
    "#Detect neighboring domains\n",
    "nbr_domians=spg.find_neighbor_clusters(target_cluster=target,\n",
    "                                   cell_id=raw.obs.index.tolist(), \n",
    "                                   x=raw.obs[\"x_array\"].tolist(), \n",
    "                                   y=raw.obs[\"y_array\"].tolist(), \n",
    "                                   pred=raw.obs[\"pred\"].tolist(),\n",
    "                                   radius=r,\n",
    "                                   ratio=1/2)\n",
    "\n",
    "nbr_domians=nbr_domians[0:3]\n",
    "de_genes_info=spg.rank_genes_groups(input_adata=raw,\n",
    "                                target_cluster=target,\n",
    "                                nbr_list=nbr_domians, \n",
    "                                label_col=\"pred\", \n",
    "                                adj_nbr=True, \n",
    "                                log=True)\n",
    "#Filter genes\n",
    "de_genes_info=de_genes_info[(de_genes_info[\"pvals_adj\"]<0.05)]\n",
    "filtered_info=de_genes_info\n",
    "filtered_info=filtered_info[(filtered_info[\"pvals_adj\"]<0.05) &\n",
    "                            (filtered_info[\"in_out_group_ratio\"]>min_in_out_group_ratio) &\n",
    "                            (filtered_info[\"in_group_fraction\"]>min_in_group_fraction) &\n",
    "                            (filtered_info[\"fold_change\"]>min_fold_change)]\n",
    "filtered_info=filtered_info.sort_values(by=\"in_group_fraction\", ascending=False)\n",
    "filtered_info[\"target_dmain\"]=target\n",
    "filtered_info[\"neighbors\"]=str(nbr_domians)\n",
    "print(\"SVGs for domain \", str(target),\":\", filtered_info[\"genes\"].tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "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>genes</th>\n",
       "      <th>in_group_fraction</th>\n",
       "      <th>out_group_fraction</th>\n",
       "      <th>in_out_group_ratio</th>\n",
       "      <th>in_group_mean_exp</th>\n",
       "      <th>out_group_mean_exp</th>\n",
       "      <th>fold_change</th>\n",
       "      <th>pvals_adj</th>\n",
       "      <th>target_dmain</th>\n",
       "      <th>neighbors</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>CAMK2N1</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.944964</td>\n",
       "      <td>1.058242</td>\n",
       "      <td>2.333675</td>\n",
       "      <td>1.578288</td>\n",
       "      <td>2.128434</td>\n",
       "      <td>1.656040e-11</td>\n",
       "      <td>0</td>\n",
       "      <td>[3, 2]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>ENC1</td>\n",
       "      <td>0.998638</td>\n",
       "      <td>0.941848</td>\n",
       "      <td>1.060295</td>\n",
       "      <td>2.457791</td>\n",
       "      <td>1.696083</td>\n",
       "      <td>2.141931</td>\n",
       "      <td>1.552131e-03</td>\n",
       "      <td>0</td>\n",
       "      <td>[3, 2]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>GPM6A</td>\n",
       "      <td>0.997275</td>\n",
       "      <td>0.922118</td>\n",
       "      <td>1.081505</td>\n",
       "      <td>2.224006</td>\n",
       "      <td>1.561187</td>\n",
       "      <td>1.940255</td>\n",
       "      <td>8.602227e-03</td>\n",
       "      <td>0</td>\n",
       "      <td>[3, 2]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>ARPP19</td>\n",
       "      <td>0.982289</td>\n",
       "      <td>0.853583</td>\n",
       "      <td>1.150784</td>\n",
       "      <td>1.889256</td>\n",
       "      <td>1.272106</td>\n",
       "      <td>1.853637</td>\n",
       "      <td>4.823349e-02</td>\n",
       "      <td>0</td>\n",
       "      <td>[3, 2]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HPCAL1</td>\n",
       "      <td>0.851499</td>\n",
       "      <td>0.465213</td>\n",
       "      <td>1.830342</td>\n",
       "      <td>1.141321</td>\n",
       "      <td>0.406338</td>\n",
       "      <td>2.085448</td>\n",
       "      <td>9.706465e-05</td>\n",
       "      <td>0</td>\n",
       "      <td>[3, 2]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     genes  in_group_fraction  out_group_fraction  in_out_group_ratio  \\\n",
       "0  CAMK2N1           1.000000            0.944964            1.058242   \n",
       "2     ENC1           0.998638            0.941848            1.060295   \n",
       "4    GPM6A           0.997275            0.922118            1.081505   \n",
       "6   ARPP19           0.982289            0.853583            1.150784   \n",
       "1   HPCAL1           0.851499            0.465213            1.830342   \n",
       "\n",
       "   in_group_mean_exp  out_group_mean_exp  fold_change     pvals_adj  \\\n",
       "0           2.333675            1.578288     2.128434  1.656040e-11   \n",
       "2           2.457791            1.696083     2.141931  1.552131e-03   \n",
       "4           2.224006            1.561187     1.940255  8.602227e-03   \n",
       "6           1.889256            1.272106     1.853637  4.823349e-02   \n",
       "1           1.141321            0.406338     2.085448  9.706465e-05   \n",
       "\n",
       "   target_dmain neighbors  \n",
       "0             0    [3, 2]  \n",
       "2             0    [3, 2]  \n",
       "4             0    [3, 2]  \n",
       "6             0    [3, 2]  \n",
       "1             0    [3, 2]  "
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filtered_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot refinedspatial domains\n",
    "color_self = clr.LinearSegmentedColormap.from_list('pink_green', ['#3AB370',\"#EAE7CC\",\"#FD1593\"], N=256)\n",
    "for g in filtered_info[\"genes\"].tolist():\n",
    "    raw.obs[\"exp\"]=raw.X[:,raw.var.index==g]\n",
    "    ax=sc.pl.scatter(raw,alpha=1,x=\"y_pixel\",y=\"x_pixel\",color=\"exp\",title=g,color_map=color_self,show=False,size=100000/raw.shape[0])\n",
    "    ax.set_aspect('equal', 'box')\n",
    "    ax.axes.invert_yaxis()\n",
    "    plt.savefig(\"./sample_results/\"+g+\".png\", dpi=600)\n",
    "    plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**CAMK2N1**![](./sample_results/CAMK2N1.png) **ENC1**![](./sample_results/ENC1.png) **GPM6A**![](./sample_results/GPM6A.png) **ARPP19**![](./sample_results/ARPP19.png) **HPCAL1**![](./sample_results/HPCAL1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Identify Meta Gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Trying to set attribute `.obs` of view, copying.\n",
      "Trying to set attribute `.obs` of view, copying.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Add gene:  MGP\n",
      "Minus gene:  FTH1\n",
      "Absolute mean change: 0.8913243\n",
      "Number of non-target spots reduced to: 1888\n",
      "===========================================================================\n",
      "Meta gene is:  GFAP+MGP-FTH1\n",
      "===========================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Trying to set attribute `.obs` of view, copying.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Add gene:  MYL9\n",
      "Minus gene:  MBP\n",
      "Absolute mean change: 2.175557\n",
      "Number of non-target spots reduced to: 563\n",
      "===========================================================================\n",
      "Meta gene is:  GFAP+MGP-FTH1+MYL9-MBP\n",
      "===========================================================================\n",
      "Add gene:  KRT8\n",
      "Minus gene:  MT-ATP6\n",
      "Absolute mean change: 2.8935516\n",
      "Number of non-target spots reduced to: 111\n",
      "===========================================================================\n",
      "Meta gene is:  GFAP+MGP-FTH1+MYL9-MBP+KRT8-MT-ATP6\n",
      "===========================================================================\n"
     ]
    }
   ],
   "source": [
    "#Use domain 2 as an example\n",
    "target=2\n",
    "meta_name, meta_exp=spg.find_meta_gene(input_adata=raw,\n",
    "                    pred=raw.obs[\"pred\"].tolist(),\n",
    "                    target_domain=target,\n",
    "                    start_gene=\"GFAP\",\n",
    "                    mean_diff=0,\n",
    "                    early_stop=True,\n",
    "                    max_iter=3,\n",
    "                    use_raw=False)\n",
    "\n",
    "raw.obs[\"meta\"]=meta_exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot meta gene\n",
    "g=\"GFAP\"\n",
    "raw.obs[\"exp\"]=raw.X[:,raw.var.index==g]\n",
    "ax=sc.pl.scatter(raw,alpha=1,x=\"y_pixel\",y=\"x_pixel\",color=\"exp\",title=g,color_map=color_self,show=False,size=100000/raw.shape[0])\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.axes.invert_yaxis()\n",
    "plt.savefig(\"./sample_results/\"+g+\".png\", dpi=600)\n",
    "plt.close()\n",
    "\n",
    "raw.obs[\"exp\"]=raw.obs[\"meta\"]\n",
    "ax=sc.pl.scatter(raw,alpha=1,x=\"y_pixel\",y=\"x_pixel\",color=\"exp\",title=meta_name,color_map=color_self,show=False,size=100000/raw.shape[0])\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.axes.invert_yaxis()\n",
    "plt.savefig(\"./sample_results/meta_gene.png\", dpi=600)\n",
    "plt.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**start**![](./sample_results/GFAP.png) **meta gene**![](./sample_results/meta_gene.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Multiple tissue sections analysis\n",
    "In this section, we show an example on how to analysis multiple adjacent tissue sections using SpaGCN."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Mouse brain anterior**![](./sample_results/MA1_lowres.png) **Mouse brain posterior**![](./sample_results/MP1_lowres.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.1 Read in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n",
      "Variable names are not unique. To make them unique, call `.var_names_make_unique`.\n"
     ]
    }
   ],
   "source": [
    "adata1=sc.read(\"./data/Mouse_brain/MA1.h5ad\")\n",
    "adata2=sc.read(\"./data/Mouse_brain/MP1.h5ad\")\n",
    "img1=cv2.imread(\"./data/Mouse_brain/MA1_histology.tif\")\n",
    "img2=cv2.imread(\"./data/Mouse_brain/MP1_histology.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.2 Extract color intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "b=49\n",
    "s=1\n",
    "x_pixel1=adata1.obs[\"x4\"].tolist()\n",
    "y_pixel1=adata1.obs[\"x5\"].tolist()\n",
    "adata1.obs[\"color\"]=spg.extract_color(x_pixel=x_pixel1, y_pixel=y_pixel1, image=img1, beta=b)\n",
    "z_scale=np.max([np.std(x_pixel1), np.std(y_pixel1)])*s\n",
    "adata1.obs[\"z\"]=(adata1.obs[\"color\"]-np.mean(adata1.obs[\"color\"]))/np.std(adata1.obs[\"color\"])*z_scale\n",
    "\n",
    "x_pixel2=adata2.obs[\"x4\"].tolist()\n",
    "y_pixel2=adata2.obs[\"x5\"].tolist()\n",
    "adata2.obs[\"color\"]=spg.extract_color(x_pixel=x_pixel2, y_pixel=y_pixel2, image=img2, beta=b)\n",
    "z_scale=np.max([np.std(x_pixel2), np.std(y_pixel2)])*s\n",
    "adata2.obs[\"z\"]=(adata2.obs[\"color\"]-np.mean(adata2.obs[\"color\"]))/np.std(adata2.obs[\"color\"])*z_scale\n",
    "del img1, img2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.3 Modify coordinates to combine 2 sections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "from anndata import AnnData\n",
    "adata1.obs[\"x_pixel\"]=x_pixel1\n",
    "adata1.obs[\"y_pixel\"]=y_pixel1\n",
    "adata2.obs[\"x_pixel\"]=x_pixel2-np.min(x_pixel2)+np.min(x_pixel1)\n",
    "adata2.obs[\"y_pixel\"]=y_pixel2-np.min(y_pixel2)+np.max(y_pixel1)\n",
    "adata1.var_names_make_unique()\n",
    "adata2.var_names_make_unique()\n",
    "adata_all=AnnData.concatenate(adata1, adata2,join='inner',batch_key=\"dataset_batch\",batch_categories=[\"0\",\"1\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.4 Integrate gene expression and histology into a Graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=np.array([adata_all.obs[\"x_pixel\"], adata_all.obs[\"y_pixel\"], adata_all.obs[\"z\"]]).T.astype(np.float32)\n",
    "adj=spg.pairwise_distance(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 8.5 Spatial domain detection using SpaGCN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Run 1: l [0.01, 1000], p [0.0, 144.17116893743565]\n",
      "Run 2: l [0.01, 500.005], p [0.0, 24.78992462158203]\n",
      "Run 3: l [0.01, 250.0075], p [0.0, 3.649960994720459]\n",
      "Run 4: l [125.00874999999999, 250.0075], p [0.4487175941467285, 3.649960994720459]\n",
      "Run 5: l [125.00874999999999, 187.508125], p [0.4487175941467285, 1.5741894245147705]\n",
      "Run 6: l [125.00874999999999, 156.2584375], p [0.4487175941467285, 0.9070142507553101]\n",
      "Run 7: l [125.00874999999999, 140.63359375], p [0.4487175941467285, 0.6537595987319946]\n",
      "Run 8: l [125.00874999999999, 132.821171875], p [0.4487175941467285, 0.5454769134521484]\n",
      "recommended l =  128.91496093749998\n"
     ]
    }
   ],
   "source": [
    "sc.pp.normalize_per_cell(adata_all, min_counts=0)\n",
    "sc.pp.log1p(adata_all)\n",
    "p=0.5 \n",
    "#Find the l value given p\n",
    "l=spg.search_l(p, adj, start=0.01, end=1000, tol=0.01, max_run=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing cluster centers with louvain, resolution =  1.0\n",
      "Epoch  0\n",
      "Epoch  10\n",
      "Epoch  20\n",
      "Epoch  30\n",
      "delta_label  0.004794180856339891 < tol  0.005\n",
      "Reach tolerance threshold. Stopping training.\n",
      "Total epoch: 37\n"
     ]
    }
   ],
   "source": [
    "res=1.0\n",
    "seed=100\n",
    "random.seed(seed)\n",
    "torch.manual_seed(seed)\n",
    "np.random.seed(seed)\n",
    "clf=spg.SpaGCN()\n",
    "clf.set_l(l)\n",
    "clf.train(adata_all,adj,init_spa=True,init=\"louvain\",res=res, tol=5e-3, lr=0.05, max_epochs=200)\n",
    "y_pred, prob=clf.predict()\n",
    "adata_all.obs[\"pred\"]= y_pred\n",
    "adata_all.obs[\"pred\"]=adata_all.obs[\"pred\"].astype('category')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors_use=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#bcbd22', '#17becf', '#aec7e8', '#ffbb78', '#98df8a', '#ff9896', '#bec1d4', '#bb7784', '#0000ff', '#111010', '#FFFF00',   '#1f77b4', '#800080', '#959595', \n",
    " '#7d87b9', '#bec1d4', '#d6bcc0', '#bb7784', '#8e063b', '#4a6fe3', '#8595e1', '#b5bbe3', '#e6afb9', '#e07b91', '#d33f6a', '#11c638', '#8dd593', '#c6dec7', '#ead3c6', '#f0b98d', '#ef9708', '#0fcfc0', '#9cded6', '#d5eae7', '#f3e1eb', '#f6c4e1', '#f79cd4']\n",
    "num_celltype=len(adata_all.obs[\"pred\"].unique())\n",
    "adata_all.uns[\"pred_colors\"]=list(colors_use[:num_celltype])\n",
    "ax=sc.pl.scatter(adata_all,alpha=1,x=\"y_pixel\",y=\"x_pixel\",color=\"pred\",show=False,size=150000/adata_all.shape[0])\n",
    "ax.set_aspect('equal', 'box')\n",
    "ax.axes.invert_yaxis()\n",
    "plt.savefig(\"./sample_results/mouse_barin_muti_sections_domains.png\", dpi=600)\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**SpaGCN mouse brain combines**![](./sample_results/mouse_barin_muti_sections_domains.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}