tutorial/tutorial_ez_mode.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1><center>SpaGCN Tutorial Easy Mode</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. Spatial domain detection using SpaGCN\n",
    "5. Identify SVGs\n",
    "6. Identify Meta Gene"
   ]
  },
  {
   "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": 5,
   "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": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.2.0'"
      ]
     },
     "execution_count": 6,
     "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": 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": [
    "\"\"\"\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",
    "#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. Spatial domain detection using SpaGCN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- s: This 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",
    "- b: This parameter determines the area of each spot when extracting color intensity.\n",
    "- 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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set coordinates\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",
    "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",
    "#Run SpaGCN\n",
    "adata.obs[\"pred\"]= spg.detect_spatial_domains_ez_mode(adata, img, x_array, y_array, x_pixel, y_pixel, n_clusters=7, histology=True, s=1, b=49, p=0.5, r_seed=100, t_seed=100, n_seed=100)\n",
    "adata.obs[\"pred\"]=adata.obs[\"pred\"].astype('category')\n",
    "#Refine domains (optional)\n",
    "#shape=\"hexagon\" for Visium data, \"square\" for ST data.\n",
    "adata.obs[\"refined_pred\"]=spg.spatial_domains_refinement_ez_mode(sample_id=adata.obs.index.tolist(), pred=adata.obs[\"pred\"].tolist(), x_array=x_array, y_array=y_array, shape=\"hexagon\")\n",
    "adata.obs[\"refined_pred\"]=adata.obs[\"refined_pred\"].astype('category')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot spatial domains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_color=[\"#F56867\",\"#FEB915\",\"#C798EE\",\"#59BE86\",\"#7495D3\",\"#D1D1D1\",\"#6D1A9C\",\"#15821E\",\"#3A84E6\",\"#997273\",\"#787878\",\"#DB4C6C\",\"#9E7A7A\",\"#554236\",\"#AF5F3C\",\"#93796C\",\"#F9BD3F\",\"#DAB370\",\"#877F6C\",\"#268785\"]\n",
    "ax=spg.plot_spatial_domains_ez_mode(adata, domain_name=\"pred\", x_name=\"y_pixel\", y_name=\"x_pixel\", plot_color=plot_color,size=150000/adata.shape[0], show=False, save=True,save_dir=\"./sample_results/pred.png\")\n",
    "ax=spg.plot_spatial_domains_ez_mode(adata, domain_name=\"refined_pred\", x_name=\"y_pixel\", y_name=\"x_pixel\", plot_color=plot_color,size=150000/adata.shape[0], show=False, save=True,save_dir=\"./sample_results/refined_pred.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Spatial Domains**![](./sample_results/pred.png) **Refined Spatial Domains**![](./sample_results/refined_pred.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Identify SVGs"
   ]
  },
  {
   "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.target: Target domain to identify SVGs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target=0\n",
    "#Read in raw data\n",
    "adata=sc.read(\"../tutorial/sample_results/results.h5ad\")\n",
    "raw=sc.read(\"../tutorial/data/151673/sample_data.h5ad\")\n",
    "#Preprocessing\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",
    "raw.X=(raw.X.A if issparse(raw.X) else raw.X)\n",
    "raw.raw=raw\n",
    "sc.pp.log1p(raw)\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",
    "filtered_info=spg.detect_SVGs_ez_mode(raw, target=0, x_name=\"x_array\", y_name=\"y_array\", domain_name=\"pred\", min_in_group_fraction=min_in_group_fraction, min_in_out_group_ratio=min_in_out_group_ratio, min_fold_change=min_fold_change)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "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": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "filtered_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "color_self = clr.LinearSegmentedColormap.from_list('pink_green', ['#3AB370',\"#EAE7CC\",\"#FD1593\"], N=256)\n",
    "spg.plot_SVGs_ez_mode(adata=raw, gene_list=filtered_info[\"genes\"].tolist(), x_name=\"y_pixel\", y_name=\"x_pixel\", plot_color=color_self, size=50, show=False, save=True, save_dir=\"./sample_results/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**CAMK2N1**![](./sample_results/CAMK2N1.png) **ENC1**![](./sample_results/ENC1.png) **GPM6A**![](./sample_results/GPM6A.png) **ARPP19**![](c/ARPP19.png) **HPCAL1**![](./sample_results/HPCAL1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Identify Meta Gene"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- target: Target domain to detect meta gene.\n",
    "- start_gene: Gene to satrt with. Candidate start genes can be found by lowering the filtering criterials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target=2\n",
    "raw.obs[\"meta\"]=spg.detect_meta_genes_ez_mode(raw, target, x_name=\"x_array\", y_name=\"y_array\", domain_name=\"pred\", start_gene=\"GFAP\", use_raw=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "spg.plot_meta_genes_ez_mode(raw, x_name=\"y_pixel\", y_name=\"x_pixel\", meta_name=\"meta\", plot_color=color_self, size=50, show=False, save=True, save_dir=\"./\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**start**![](./sample_results/GFAP.png) **meta gene**![](./sample_results/meta_gene.png)"
   ]
  }
 ],
 "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
}