{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cell-cell connection estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimation aims to process on spatially resolved ligand-receptors based on SPC and calculate cell cell interaction strength."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We modified [NICHES](https://github.com/msraredon/NICHES) for the 3D data calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First import all necessary library and custom modified NICHES_3D.R library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(Seurat)\n",
    "library(SeuratObject)\n",
    "library(SeuratDisk)\n",
    "library(dplyr)\n",
    "library(tidyr)\n",
    "library(reshape2)\n",
    "library(NICHES)\n",
    "library(stringr)\n",
    "\n",
    "source(\"./NICHES_3D.R\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in RDS formated Seurat object which include **x**, **y**, **z** 3D coordinates in obj$meta.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj <- readRDS('WT.rds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To prepare the ligand-receptor pairs for planarian, we used blast hit to [CellChat human ligand-receptor pairs](https://github.com/jinworks/CellChat)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "custom_lrdb <- read.table('./CellChatDB.human.interaction.planarian.txt', sep='\\t', header=T)\n",
    "custom_lrdb <- custom_lrdb[c('gene_id.ligand', 'gene_id.receptor')]\n",
    "custom_lrdb <- na.omit(custom_lrdb)\n",
    "custom_lrdb <- custom_lrdb[!duplicated(custom_lrdb), ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then run the NICHES pipeline with spatial cell to cell mode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results <- RunNICHES(obj,\n",
    "                     assay='SCT',\n",
    "                     LR.database=\"custom\",\n",
    "                     min.cells.per.ident = 3,\n",
    "                     min.cells.per.gene = 3,\n",
    "                     meta.data.to.map = colnames(obj@meta.data),\n",
    "                     position.x = 'x',\n",
    "                     position.y = 'y',\n",
    "                     position.z = 'z',\n",
    "                     cell_types = 'spc_cluster', \n",
    "                     custom_LR_database = custom_lrdb, \n",
    "                     k = NULL,\n",
    "                     rad.set = 60, # median_mindis * 5, five layer cells\n",
    "                     blend = 'mean',\n",
    "                     CellToCell = F,\n",
    "                     CellToSystem = F,\n",
    "                     SystemToCell = F,\n",
    "                     CellToCellSpatial = T,\n",
    "                     CellToNeighborhood = F,\n",
    "                     NeighborhoodToCell = F,\n",
    "                     output_format = \"seurat\"\n",
    "                     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj.result <- results$CellToCellSpatial\n",
    "obj.loom <- as.loom(obj.result, filename = 'WT_niches.loom', overwrite = T)\n",
    "obj.loom$close_all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will generate a matrix in loom format with lr pairs as feature and cell pairs as obs. With the calculated interaction strength and cell location, connections among cells were constucted with cells as nodes and sender-reseiver relationships as edges. \n",
    "\n",
    "The direction of the connection was calculated as the angle of each edge in the 3D space, which was reflected by the orientation of the edge. The entropy of the connection was estimated based on the degree of randomness of the connection orientations. The entropy was used to infer the region in which highly synergistic cell signaling activities occurred during the remodeling of regeneration following: \n",
    " \n",
    "$$\n",
    "Pr⁡[X=s] ≅ freq(s) = \\sum_{X=s} X_i \\div \\sum_{all directions} X_i\n",
    "$$\n",
    "$$\n",
    "estimate(H(x)) = - \\sum_{all directions} freq(s) * log⁡(freq(s)) .\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}