Quenching_Single_Cell.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single Cell Analysis of Quenched Slides\n",
    "\n",
    "**Question:** How does autofluorescence change over repeated quenching reactions?\n",
    "\n",
    "**Samples:** Three **normal pancreas** tissue sections were processed with the cyclif IF protocol, then repeatedly quenched for 15 minutes and imaged. No anitbody was applied. Quenching solution contained 3%, 4.5% and 6% H2O2, but no difference in quenching speed or completeness was observed between the conditions.\n",
    "\n",
    "Three adjacent sections from **Biomax FDA808l-2** :  Multiple tumor (24 organs )and normal (6 organs) tissue array, with grade and TNM/Stage,72 cases/72 cores(core size 1.5mm), replacing FDA808k-2, were quenched and imaged repeatedly. (10 rounds)\n",
    "\n",
    "HER2A TMA: **Biomax Her2a** : Her-2 IHC Bio-marker control TMA, with breast invasive ductal carcinoma, including TNM, clinical stage and pathology grade, 10 cases/ 10 cores, replaced by Her2b. Scenes #1, 4,5 8,9 good scenes til R10  #10, 7, 6 bad after R6  #2, 3 missing after R6; 11 missing after R6\n",
    "\n",
    "**Method**: Slides were imaged as described above, registered and segmented based on DAPI stain. Nuclei were inflated by 5 pixels to roughly capture cytoplasmic signal.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib\n",
    "import seaborn as sns\n",
    "matplotlib.rcParams.update({'font.size': 15})\n",
    "import umap\n",
    "from sklearn.preprocessing import scale\n",
    "import scanpy as sc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set location of files\n",
    "os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cycIF_ValidationStudies/cycIF_Validation')\n",
    "codedir = os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import and Filter Data\n",
    "\n",
    "## Normal Tissue TMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "ls_sample = ['808L2-053','808L2-054','808L2-055'] #\n",
    "#ls_sample = ['HER2A-60','HER2A-61','HER2A-62']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi = pd.DataFrame()\n",
    "df_xy = pd.DataFrame()\n",
    "for s_sample in ls_sample:\n",
    "    df_mi = df_mi.append(pd.read_csv(f'{codedir}/Data/Quench/features_{s_sample}_patched_MeanIntensity_Shape__registeredimages.csv',index_col=0))\n",
    "    df_xy = df_xy.append(pd.read_csv(f'{codedir}/Data/Quench/features_{s_sample}_CentroidXY.csv',index_col=0))\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filtering out edge cells from adjacent TMA cores\n",
    "ls_drop = ['808L2-053_scene066',#wrong name\n",
    "          ##'808L2-055_scene023', '808L2-055_scene052','808L2-055_scene056', '808L2-055_scene066', #missing mean int\n",
    "           #'808L2-055_scene071' #missing df_xy\n",
    "          ] \n",
    "\n",
    "#### BOTTOM ############\n",
    "d_filter = {'808L2-053_scene006':(df_xy.DAPI_Y > 6000),'808L2-053_scene007':(df_xy.DAPI_Y > 7300),\n",
    "            '808L2-053_scene019':(df_xy.DAPI_Y > 7000),'808L2-053_scene028':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-053_scene020':(df_xy.DAPI_Y > 6000),'808L2-053_scene021':(df_xy.DAPI_Y > 6000),\n",
    "            '808L2-053_scene027':(df_xy.DAPI_Y > 6000),'808L2-053_scene029':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-053_scene030':(df_xy.DAPI_Y > 8000),'808L2-053_scene036':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-053_scene037':(df_xy.DAPI_Y > 8000),'808L2-053_scene038':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-053_scene046':(df_xy.DAPI_Y > 8000),'808L2-053_scene052':(df_xy.DAPI_Y > 6000),\n",
    "           '808L2-053_scene053':(df_xy.DAPI_Y > 8000),'808L2-053_scene058':(df_xy.DAPI_Y > 6000),\n",
    "           '808L2-053_scene059':(df_xy.DAPI_Y > 8000),\n",
    "           #54\n",
    "           '808L2-054_scene012':(df_xy.DAPI_Y > 8000),'808L2-054_scene013':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-054_scene019':(df_xy.DAPI_Y > 8000),'808L2-054_scene026':(df_xy.DAPI_Y > 6000),\n",
    "           '808L2-054_scene027':(df_xy.DAPI_Y > 6000),'808L2-054_scene033':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-054_scene035':(df_xy.DAPI_Y > 8000),'808L2-054_scene036':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-054_scene041':(df_xy.DAPI_Y > 8000),'808L2-054_scene043':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-054_scene049':(df_xy.DAPI_Y > 7000),'808L2-054_scene050':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-054_scene052':(df_xy.DAPI_Y > 8000),'808L2-054_scene056':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-054_scene058':(df_xy.DAPI_Y > 8000),'808L2-054_scene034':(df_xy.DAPI_Y > 8000),\n",
    "            '808L2-054_scene044':(df_xy.DAPI_Y > 8000),'808L2-054_scene057':(df_xy.DAPI_Y > 8000),\n",
    "           #55\n",
    "           '808L2-055_scene012':(df_xy.DAPI_Y > 8000),'808L2-055_scene019':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-055_scene021':(df_xy.DAPI_Y > 8000),'808L2-055_scene025':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-055_scene030':(df_xy.DAPI_Y > 8000),'808L2-055_scene031':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-055_scene035':(df_xy.DAPI_Y > 8000),'808L2-055_scene040':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-055_scene041':(df_xy.DAPI_Y > 7000),'808L2-055_scene044':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-055_scene046':(df_xy.DAPI_Y > 8000),'808L2-055_scene054':(df_xy.DAPI_Y > 8000),\n",
    "           '808L2-055_scene055':(df_xy.DAPI_Y > 8000),'808L2-055_scene059':(df_xy.DAPI_Y > 8000),\n",
    "            '808L2-055_scene037':(df_xy.DAPI_Y > 6000),\n",
    "           }\n",
    "ls_filter_all = []\n",
    "for s_scene, filtercon in d_filter.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "##### LEFT ############   \n",
    "d_filter2 = {'808L2-053_scene011':(df_xy.DAPI_X < 1500),'808L2-053_scene019':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene020':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene028':(df_xy.DAPI_X < 1500),'808L2-053_scene029':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene030':(df_xy.DAPI_X < 1500),'808L2-053_scene037':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene038':(df_xy.DAPI_X < 1500),'808L2-053_scene046':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene047':(df_xy.DAPI_X < 1500),'808L2-053_scene048':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene053':(df_xy.DAPI_X < 1500),'808L2-053_scene057':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene064':(df_xy.DAPI_X < 1500),'808L2-053_scene065':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-053_scene072':(df_xy.DAPI_X < 1500),'808L2-053_scene056':(df_xy.DAPI_X < 1500),\n",
    "             '808L2-053_scene069':(df_xy.DAPI_X < 1500),\n",
    "            #54\n",
    "             '808L2-054_scene007':(df_xy.DAPI_X < 1500),'808L2-054_scene020':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-054_scene024':(df_xy.DAPI_X < 1500),'808L2-054_scene025':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-054_scene027':(df_xy.DAPI_X < 1500),'808L2-054_scene029':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-054_scene030':(df_xy.DAPI_X < 900),'808L2-054_scene038':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-054_scene047':(df_xy.DAPI_X < 1500),'808L2-054_scene048':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-054_scene049':(df_xy.DAPI_X < 1500),'808L2-054_scene058':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-054_scene068':(df_xy.DAPI_X < 1500),'808L2-054_scene069':(df_xy.DAPI_X < 1500),\n",
    "            #55\n",
    "            '808L2-055_scene003':(df_xy.DAPI_X < 1500),'808L2-055_scene004':(df_xy.DAPI_X < 1500),\n",
    "            '808L2-055_scene020':(df_xy.DAPI_X < 1000),'808L2-055_scene024':(df_xy.DAPI_X < 1000),\n",
    "            '808L2-055_scene026':(df_xy.DAPI_X < 1000),'808L2-055_scene030':(df_xy.DAPI_X < 1000),\n",
    "            '808L2-055_scene031':(df_xy.DAPI_X < 1000),'808L2-055_scene032':(df_xy.DAPI_X < 1000),\n",
    "            '808L2-055_scene053':(df_xy.DAPI_X < 1000),'808L2-055_scene054':(df_xy.DAPI_X < 1000),\n",
    "            '808L2-055_scene060':(df_xy.DAPI_X < 1000),'808L2-055_scene068':(df_xy.DAPI_X < 1000),\n",
    "            '808L2-055_scene069':(df_xy.DAPI_X < 1000),'808L2-055_scene070':(df_xy.DAPI_X < 1000),\n",
    "             '808L2-055_scene047':(df_xy.DAPI_X < 1000),'808L2-055_scene061':(df_xy.DAPI_X < 1000),\n",
    "            }\n",
    "for s_scene, filtercon in d_filter2.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "##### RIGHT ############\n",
    "d_filter3 = {'808L2-053_scene001':(df_xy.DAPI_X > 8000),'808L2-053_scene002':(df_xy.DAPI_X > 8000),\n",
    "             '808L2-053_scene019':(df_xy.DAPI_X > 8000),'808L2-053_scene027':(df_xy.DAPI_X > 6500),\n",
    "            '808L2-053_scene028':(df_xy.DAPI_X > 8000),'808L2-053_scene029':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-053_scene032':(df_xy.DAPI_X > 8000),'808L2-053_scene037':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-053_scene038':(df_xy.DAPI_X > 8000),'808L2-053_scene046':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-053_scene047':(df_xy.DAPI_X > 8000),'808L2-053_scene057':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-053_scene064':(df_xy.DAPI_X > 8500),'808L2-053_scene067':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-053_scene071':(df_xy.DAPI_X > 6000),'808L2-053_scene052':(df_xy.DAPI_X > 7000),\n",
    "            '808L2-053_scene054':(df_xy.DAPI_X > 8000),'808L2-053_scene045':(df_xy.DAPI_X > 6000),\n",
    "             '808L2-053_scene056':(df_xy.DAPI_X > 8000),'808L2-053_scene063':(df_xy.DAPI_X > 8000),\n",
    "             '808L2-053_scene066':(df_xy.DAPI_X > 6000),\n",
    "            #54\n",
    "            '808L2-054_scene007':(df_xy.DAPI_X > 8000),'808L2-054_scene010':(df_xy.DAPI_X > 6000),\n",
    "            '808L2-054_scene019':(df_xy.DAPI_X > 8000),'808L2-054_scene024':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-054_scene028':(df_xy.DAPI_X > 8000),'808L2-054_scene029':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-054_scene030':(df_xy.DAPI_X > 8000),'808L2-054_scene033':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-054_scene038':(df_xy.DAPI_X > 8000),'808L2-054_scene047':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-054_scene048':(df_xy.DAPI_X > 8000),'808L2-054_scene049':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-054_scene058':(df_xy.DAPI_X > 8000),'808L2-054_scene059':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-054_scene060':(df_xy.DAPI_X > 8000),'808L2-054_scene065':(df_xy.DAPI_X > 7000),\n",
    "            #55\n",
    "            '808L2-055_scene003':(df_xy.DAPI_X > 8000),'808L2-055_scene019':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-055_scene024':(df_xy.DAPI_X > 8000),'808L2-055_scene025':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-055_scene026':(df_xy.DAPI_X > 8000),'808L2-055_scene031':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-055_scene038':(df_xy.DAPI_X > 7000),'808L2-055_scene046':(df_xy.DAPI_X > 8000),\n",
    "            '808L2-055_scene053':(df_xy.DAPI_X > 8000),'808L2-055_scene067':(df_xy.DAPI_X > 7000),\n",
    "            '808L2-055_scene068':(df_xy.DAPI_X > 8000),'808L2-055_scene033':(df_xy.DAPI_X > 8000),\n",
    "             '808L2-055_scene047':(df_xy.DAPI_X > 8000),\n",
    "            }\n",
    "for s_scene, filtercon in d_filter3.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "    \n",
    "### TOP ############\n",
    "d_filter4 = {'808L2-053_scene014':(df_xy.DAPI_Y < 2000), '808L2-053_scene028':(df_xy.DAPI_Y < 1500),\n",
    "             '808L2-053_scene029':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-053_scene036':(df_xy.DAPI_Y < 1500),'808L2-053_scene037':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-053_scene038':(df_xy.DAPI_Y < 1500),'808L2-053_scene053':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-053_scene060':(df_xy.DAPI_Y < 1500),'808L2-053_scene061':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-053_scene062':(df_xy.DAPI_Y < 1200),'808L2-053_scene067':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-053_scene068':(df_xy.DAPI_Y < 1500),'808L2-053_scene066':(df_xy.DAPI_Y < 1500),\n",
    "            #54\n",
    "            '808L2-054_scene011':(df_xy.DAPI_Y < 1500),'808L2-054_scene012':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene033':(df_xy.DAPI_Y < 1500),'808L2-054_scene035':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene036':(df_xy.DAPI_Y < 1500),'808L2-054_scene041':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene042':(df_xy.DAPI_Y < 1500),'808L2-054_scene043':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene045':(df_xy.DAPI_Y < 1500),'808L2-054_scene050':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene052':(df_xy.DAPI_Y < 1500),'808L2-054_scene056':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene058':(df_xy.DAPI_Y < 1500),'808L2-054_scene065':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-054_scene066':(df_xy.DAPI_Y < 1500),'808L2-054_scene034':(df_xy.DAPI_Y < 1500),\n",
    "             '808L2-054_scene044':(df_xy.DAPI_Y < 1500),'808L2-054_scene057':(df_xy.DAPI_Y < 1500),\n",
    "            #55\n",
    "            '808L2-055_scene010':(df_xy.DAPI_Y < 1500),'808L2-055_scene011':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene012':(df_xy.DAPI_Y < 2000),'808L2-055_scene015':(df_xy.DAPI_Y < 2000),\n",
    "            '808L2-055_scene031':(df_xy.DAPI_Y < 2000),'808L2-055_scene035':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene036':(df_xy.DAPI_Y < 1500),'808L2-055_scene044':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene046':(df_xy.DAPI_Y < 1500),'808L2-055_scene049':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene050':(df_xy.DAPI_Y < 1500),'808L2-055_scene054':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene055':(df_xy.DAPI_Y < 1500),'808L2-055_scene059':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene060':(df_xy.DAPI_Y < 1500),'808L2-055_scene063':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene064':(df_xy.DAPI_Y < 1500),'808L2-055_scene065':(df_xy.DAPI_Y < 1500),\n",
    "            '808L2-055_scene068':(df_xy.DAPI_Y < 1500),'808L2-055_scene069':(df_xy.DAPI_Y < 2000),\n",
    "            }\n",
    "for s_scene, filtercon in d_filter4.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "    \n",
    "ls_drop_filter = df_xy[df_xy.slide_scene.isin(ls_drop)].index.tolist() #[]#\n",
    "# drop bad scenes\n",
    "ls_filter_all = ls_filter_all + ls_drop_filter #ls_filter\n",
    "#filter edge\n",
    "#ls_filter_all = ls_filter_all + df_edge.index.tolist()\n",
    "df_filter_mi=df_mi[~df_mi.index.isin(ls_filter_all)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for s_scene in ls_drop:\n",
    "    #print(f'sbatch registration_matlab_{s_scene.replace(\"-\",\"\").replace(\"scene0\",\"Scene\")}.sbatch')\n",
    "    #print(f'sbatch segmentation_nuccell_{s_scene.replace(\"-\",\"\")}.sbatch')\n",
    "    #print(f'sbatch segmentation_match_{s_scene.replace(\"-\",\"\")}.sbatch')\n",
    "ls_drop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot spatially  \n",
    "\n",
    "for s_scene in sorted(set(df_xy.slide_scene)):\n",
    "    s_scene = '808L2-053_scene066'\n",
    "    #s_scene = 'HER2A-60_scene011'\n",
    "    if len(df_xy[df_xy.slide_scene==s_scene])>=1:\n",
    "        fig,ax = plt.subplots(dpi=100)\n",
    "        #plot negative cells\n",
    "        ax.scatter(data=df_xy[df_xy.slide_scene==s_scene],x='DAPI_X',y='DAPI_Y',color='silver',s=0.3,label=f'')\n",
    "        if len(df_filter_mi[(df_filter_mi.slide_scene==s_scene)])>=1:\n",
    "                #plot positive cells\n",
    "                ls_pos = df_filter_mi[(df_filter_mi.slide_scene==s_scene)].index\n",
    "                ax.scatter(data=df_xy.loc[ls_pos],x='DAPI_X',y='DAPI_Y',s=0.3,color='darkblue') #label=f'{s_color}',               \n",
    "        ax.set_title(f'{s_scene}')\n",
    "        ax.axis('equal')\n",
    "        ax.set_ylim(ax.get_ylim()[::-1])\n",
    "        #ax[idx].set_xticklabels('')\n",
    "        #ax[idx].set_yticklabels('')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### rename scenes\n",
    "\n",
    "match across TMAS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_scene = sorted(set([item.split('_')[1] for item in df_xy.slide_scene]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#non matching scenes\n",
    "d_replace = dict(zip(['808L2-053_' + item for item in ls_scene[17:65]],['808L2-053_' + item for item in ls_scene[18:66]]))\n",
    "d_replace.update(dict(zip(['808L2-055_' + item for item in ls_scene[66:73]],['808L2-055_' + item for item in ls_scene[65:72]])))\n",
    "d_replace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### drop theis scene from df_xy before renamimg\n",
    "df_xy = df_xy[~df_xy.slide_scene.isin(ls_drop)]\n",
    "df_xy['slide_scene_new'] = df_xy.slide_scene.replace(d_replace)\n",
    "df_mi['slide_scene_new'] = df_mi.slide_scene.replace(d_replace)\n",
    "#df_filter_mi['slide_scene_new'] = df_filter_mi.slide_scene.replace(d_replace)\n",
    "df_filter_mi = df_filter_mi.merge(df_mi.loc[:,'slide_scene_new'],left_index=True,right_index=True) #,inplace=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#do scenes match\n",
    "for s_scene in sorted(set([item.split('_')[1] for item in df_xy.slide_scene])):\n",
    "    fig,ax = plt.subplots(1,3,dpi=100, figsize=(12,4))\n",
    "    for idx, s_slide in enumerate(ls_sample):\n",
    "        s_scene ='scene011'\n",
    "        #plot negative cells\n",
    "        if len(df_xy[df_xy.slide_scene_new==f'{s_slide}_{s_scene}'])>=1:\n",
    "            ax[idx].scatter(data=df_xy[df_xy.slide_scene_new==f'{s_slide}_{s_scene}'],x='DAPI_X',y='DAPI_Y',color='darkblue',s=0.3,label=f'')              \n",
    "        ax[idx].set_title(f'{s_slide}_{s_scene}')\n",
    "        ax[idx].axis('equal')\n",
    "        ax[idx].set_ylim(ax[idx].get_ylim()[::-1])\n",
    "        ax[idx].set_xticklabels('')\n",
    "        ax[idx].set_yticklabels('')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fix index\n",
    "df_filter_mi.index = [df_filter_mi.loc[item,'slide_scene_new'] + '_' + item.split('_')[2] for item in df_filter_mi.index]\n",
    "df_xy.index = [df_xy.loc[item,'slide_scene_new'] + '_' + item.split('_')[2] for item in df_xy.index]\n",
    "df_filter_mi['scene'] = [item.split('_')[1] for item in df_filter_mi.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_filter_mi.loc[:,((df_filter_mi.dtypes!='float64') | (df_filter_mi.columns.isin(['nuclei_eccentricity'])))]\n",
    "df_filter_mi['scene'] = [item.split('_')[1] for item in df_filter_mi.index]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save tissue loss\n",
    "s_slides = s_sample.split('-')[0] #'HER2A'#'808L2'# \n",
    "s_out = f'{codedir}/Data/Quench/features_{s_slides}_TissueLoss.csv'\n",
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_filter_mi.loc[:,((df_filter_mi.dtypes!='float64') | (df_filter_mi.columns.isin(['nuclei_eccentricity'])))].to_csv(s_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter R10 dapi (808)\n",
    "if s_slides == '808L2':\n",
    "    df_filter_mi = df_filter_mi[df_filter_mi.DAPI10_nuclei_thresh.fillna(False)]\n",
    "# filter R6 dapi (HER2A)\n",
    "if s_slides == 'HER2A':\n",
    "    df_filter_mi = df_filter_mi[df_filter_mi.DAPI6_nuclei_thresh.fillna(False)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save for clustering\n",
    "\n",
    "for s_col in ['c2_perinuc5','c3_perinuc5','c2_nuclei','c3_nuclei']:\n",
    "    if s_slides == '808L2':\n",
    "        ls_drop_scene = ['scene006','scene012','scene013','scene016', 'scene072'] #scene6 lost in 55, 72 mis registered in 53, others few cells\n",
    "        ls_drop_cols = []\n",
    "    if s_slides == 'HER2A':\n",
    "        ls_drop_scene = ['scene011'] #HER2A\n",
    "        ls_drop_cols = ['R7c2_perinuc5', 'R8c2_perinuc5', 'R9c2_perinuc5','R10c2_perinuc5',\n",
    "                       'R7c3_perinuc5', 'R8c3_perinuc5', 'R9c3_perinuc5','R10c3_perinuc5',\n",
    "                       'R7c2_nuclei', 'R8c2_nuclei', 'R9c2_nuclei','R10c2_nuclei',\n",
    "                       'R7c3_nuclei', 'R8c3_nuclei', 'R9c3_nuclei','R10c3_nuclei'] #HER2A goot until R6\n",
    "\n",
    "    s_out = f'{codedir}/Data/Quench/features_{s_slides}_{s_col}_MeanIntensity.csv'\n",
    "    if not os.path.exists(s_out):\n",
    "        print('saving csv')\n",
    "        df_filter_mi.loc[~df_filter_mi.scene.isin(ls_drop_scene),((df_filter_mi.columns.str.contains(s_col)) & (~df_filter_mi.columns.isin(ls_drop_cols)))].dropna().to_csv(s_out)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save corrected xy\n",
    "if not os.path.exists(f'{codedir}/Data/Quench/features_{s_slides}_CentroidXY.csv'):\n",
    "    print('saving csv')\n",
    "    df_xy.to_csv(f'{codedir}/Data/Quench/features_{s_slides}_CentroidXY.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normal Pancreas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import single cell intensity data\n",
    "ls_sample = ['4165NPanc-74','4165NPanc-76','4165NPanc-78']\n",
    "df_data = pd.DataFrame()\n",
    "df_xy=pd.DataFrame()\n",
    "for s_sample in ls_sample:\n",
    "    print(s_sample)\n",
    "    df = pd.read_csv(f'{codedir}/Data/Quench/features_{s_sample}_FilteredMeanIntensity_DAPI6_DAPI1.csv',index_col=0)\n",
    "    df['sample'] = s_sample\n",
    "    df_data = df_data.append(df)\n",
    "    df_xy = df_xy.append(pd.read_csv(f'{codedir}/Data/Quench/features_{s_sample}_CentroidXY.csv',index_col=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select nucleus and 5 pixels around nucleus (perinuc5) mean intensity features\n",
    "ls_cyto = sorted(df_data.columns[(df_data.columns.str.contains('perinuc5')) & (~df_data.columns.str.contains('DAPI')) & (~df_data.columns.str.contains('R6'))]) \n",
    "ls_nuc = [item.split('_')[0] + '_nuclei' for item in ls_cyto] \n",
    "ls_cyto = ls_cyto + ['sample']\n",
    "df_cyto = df_data.loc[:,ls_cyto].dropna()\n",
    "df_nuc = df_data.loc[:,ls_nuc].dropna()\n",
    "df_both = df_data.loc[:,ls_cyto+ls_nuc].dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analyze only channel 2 nd 3 (488 and 555)\n",
    "df_af = df_both.loc[:,(df_both.columns.str.contains('c2') |  df_both.columns.str.contains('sample'))] #df_both.columns.str.contains('c3') | #shading artifact\n",
    "\n",
    "#sort column names\n",
    "df_af = df_af.loc[:,sorted(df_af.columns)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save as csv for scanpy\n",
    "df_af.columns = [item.replace('_nuclei','_nuc') for item in df_af.columns]\n",
    "df_af.columns = [item.replace('_perinuc5','_cyt') for item in df_af.columns]\n",
    "df_af.drop('sample',axis=1).to_csv(f'{codedir}/Data/Quench/features_4165NPanc_AF_scanpy.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scanpy\n",
    "\n",
    "single cell analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#generate UMAP\n",
    "s_sample = 'features_4165NPanc_AF_scanpy'\n",
    "s_sample = 'HER2A_c3_nuclei_MeanIntensity'\n",
    "s_sample = '808L2_c2_nuclei_MeanIntensity' # perinuc5 # nuclei #c2 c3\n",
    "print(s_sample)\n",
    "\n",
    "s_slide = s_sample.split('_')[0]\n",
    "s_channel = s_sample.split('_')[1]\n",
    "s_nuc = s_sample.split('_')[2]\n",
    "filename = f'{codedir}/Data/Quench/features_{s_sample}.csv'\n",
    "print(filename)\n",
    "adata = sc.read_csv(filename)\n",
    "\n",
    "#adata.obs\n",
    "if s_slide == '4165NPanc':\n",
    "    adata.obs['tissue'] = [item.split('_scene')[0].split('4165')[1] for item in adata.obs.index]\n",
    "else:\n",
    "    adata.obs['tissue'] = [item.split('_scene')[0] for item in adata.obs.index]\n",
    "#log transform, batch correct reduce dimensionality (PCA)\n",
    "sc.pp.log1p(adata)\n",
    "adata.raw = adata\n",
    "\n",
    "#remove batch effect\n",
    "sc.pp.combat(adata,key='tissue')\n",
    "\n",
    "#reduce dimensionality\n",
    "sc.tl.pca(adata, svd_solver='auto')\n",
    "if s_slide == '808L2':\n",
    "    n_neighbors = 25 #25 for 808 #20 for HER2\n",
    "elif s_slide == 'HER2A':\n",
    "    n_neighbors = 20 \n",
    "else:\n",
    "    n_neighbors=10\n",
    "print(f'Finding {n_neighbors} neighbors')\n",
    "sc.pp.neighbors(adata, n_neighbors=n_neighbors) #10 for 1 tissue\n",
    "print('calc Umap')\n",
    "sc.tl.umap(adata)\n",
    "results_file = f'{codedir}/Data/Quench/{s_slide}_{s_channel}_{s_nuc}_norm_{n_neighbors}-5s.h5ad'  # the file that will store the analysis results\n",
    "adata.write(results_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#umap run_scanpy.py\n",
    "'''\n",
    "sc.pp.neighbors(adata, n_neighbors=25) #10 for 1 tissue (4165), 25 for TMA\n",
    "sc.tl.umap(adata)\n",
    "results_file = f'{codedir}/Data/{s_slide}_AF_norm.h5ad'  # the file that will store the analysis results\n",
    "#adata.write(results_file)\n",
    "#results_file = f'{codedir}/Data/{s_slide}_AF_raw.h5ad'  # the file that will store the analysis results\n",
    "#rdata.write(results_file)\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# load umap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = 'features_4165NPanc_AF_scanpy'\n",
    "s_sample = '808L2_c2_nuclei_MeanIntensity'\n",
    "#s_sample = 'HER2A_c3_nuclei_MeanIntensity'\n",
    "s_slide = s_sample.split('_')[0]\n",
    "s_sample = s_sample.split('_MeanIntensity')[0]\n",
    "n_neighbors= 25 #10 #25 for 808 #20 for HER2A\n",
    "if s_slide == '808L2':\n",
    "    if n_neighbors==25:\n",
    "        filename = f'{codedir}/Data/{s_sample}_norm_{n_neighbors}-5s.h5ad'\n",
    "    else:\n",
    "        filename = f'{codedir}/Data/{s_slide}_AF_norm_{n_neighbors}-5s.h5ad'\n",
    "else:\n",
    "    filename = f'{codedir}/Data/{s_sample}_norm_{n_neighbors}-5s.h5ad'\n",
    "print(filename)\n",
    "rawfilename = filename\n",
    "adata = sc.read_h5ad(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#umap plot\n",
    "figname = f'{codedir}/Figures/{s_sample}_UmapBatch_{n_neighbors}.png'\n",
    "fig, ax = plt.subplots(figsize = (3,3),dpi=300)\n",
    "sc.pl.umap(adata, color='tissue',ax=ax,title='Normalized')\n",
    "plt.show(fig)\n",
    "plt.tight_layout()\n",
    "#fig.savefig(figname, dpi=150)\n",
    "fig.savefig(figname,bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# umap plus expression\n",
    "sc.set_figure_params(fontsize=16)\n",
    "for s_color in sorted(adata.var.index):\n",
    "    #s_color = 'R0c2_nuclei'#'R2c2_nuc'\n",
    "    figname = f'{codedir}/Figures/{s_sample}_{n_neighbors}_{s_color}_Umap.png'\n",
    "    fig, ax = plt.subplots(figsize = (4,4),dpi=300)\n",
    "    sc.pl.umap(adata, color=s_color, wspace=0.2,title=s_color,ax=ax,vmin=6,vmax=8) #size = 1.5, vmin='p1.5',vmax='p98.5'\n",
    "    plt.tight_layout()\n",
    "    plt.show(fig)\n",
    "    fig.savefig(figname,bbox_inches='tight')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load annotation\n",
    "df_a = pd.read_csv(f'{codedir}/Data/FDA808l-2 specs.csv')\n",
    "#grade\n",
    "df_a['Grade'] = df_a.Grade.replace('-',np.nan)\n",
    "#stage\n",
    "d_replace = {'-':np.nan, 'IA':'I', 'III':\"III\", 'IIA':\"II\", 'IIB':'II', 'IIIA':\"III\", 'IB':\"I\",\n",
    "       ' IIB G2':\"II\", 'IIIB':\"III\", 'IVB':'IV', 'IIA G3':\"II\", 'IIB G3 ':\"II\", 'IIB G2':\"II\"}\n",
    "df_a['Stage'] = df_a.Stage.replace(d_replace)\n",
    "df_a.loc[:,'Pathology'] = df_a.loc[:,'Pathology diagnosis']\n",
    "df_a.loc[:,'Age_Q'] = pd.qcut(df_a.Age,4,labels=[\"Q1\", \"Q2\", \"Q3\",\"Q4\"])\n",
    "# add scene\n",
    "adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_type in ['Type','Grade','Stage','Sex','Organ','Pathology','Age_Q']: #\n",
    "    figname = f'{codedir}/Figures/{s_sample}_{n_neighbors}_{s_type}_Umap.png'\n",
    "    fig, ax = plt.subplots(figsize = (5,4),dpi=150)\n",
    "    adata.obs[s_type] = adata.obs['scene'].map(dict(zip(df_a.Scene,df_a.loc[:,s_type])))\n",
    "    sc.pl.umap(adata, color=s_type, wspace=0.2,title=s_type,ax=ax) \n",
    "    plt.tight_layout()\n",
    "    plt.show(fig)\n",
    "    fig.savefig(figname,bbox_inches='tight')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# umap plus annotation\n",
    "'''\n",
    "s_color = 'R10c2_nuclei'#'R5c2_nuclei'#'R5c2_nuc'\n",
    "sc.set_figure_params(fontsize=16)\n",
    "figname = f'{codedir}/Figures/{s_slide}_{s_color}_Umap{n_neighbors}.png'\n",
    "fig, ax = plt.subplots(figsize = (5,4),dpi=150)\n",
    "sc.pl.umap(adata, color=s_color, wspace=0.2,title='AF488: 60 min',ax=ax,vmin=5,vmax=9)\n",
    "plt.tight_layout()\n",
    "plt.show(fig)\n",
    "fig.savefig(figname,bbox_inches='tight')\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cluster\n",
    "if s_slide == '4165NPanc':\n",
    "    resolution = 0.3\n",
    "else:\n",
    "    resolution = 0.2\n",
    "print(resolution)\n",
    "if adata.obs.columns.isin(['leiden']).sum() == 0:\n",
    "    sc.tl.leiden(adata, resolution = resolution)\n",
    "# umap plus clustering result \n",
    "sc.set_figure_params(fontsize=16)\n",
    "figname = f'{codedir}/Figures/{s_sample}_UmapLeiden_neigh{n_neighbors}.png'\n",
    "fig, ax = plt.subplots(figsize = (4,4),dpi=150)\n",
    "sc.pl.umap(adata, color='leiden',wspace=0.2,title='Clustering w/ Leiden algorithm',ax=ax) #size=1,\n",
    "plt.tight_layout()\n",
    "plt.show(fig)\n",
    "fig.savefig(figname,bbox_inches='tight')\n",
    "#save results\n",
    "#resolution = 0.2\n",
    "results_file = f'{codedir}/Data/{s_sample}_norm_{n_neighbors}_leiden{resolution}.h5ad'  # the file that will store the analysis results\n",
    "if not os.path.exists(results_file):\n",
    "    adata.write(results_file)\n",
    "#save as csv\n",
    "#adata = sc.read_h5ad(results_file)\n",
    "df_data = pd.DataFrame(adata.X, columns = adata.var.index, index = adata.obs.index)\n",
    "#df_raw = pd.DataFrame(rdata.raw.X, columns = adata.var.index, index = adata.obs.index)\n",
    "#save to csv\n",
    "df_data.merge(adata.obs,left_index=True,right_index=True).to_csv(f'{codedir}/Data/{s_sample}_norm_{n_neighbors}_leiden{resolution}.csv')\n",
    "#df_raw.merge(adata.obs,left_index=True,right_index=True).to_csv(f'{codedir}/Data/{s_slide}_AF_raw_{n_neighbors}_leiden{resolution}.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '808L2_c2_nuclei'# = 'HER2A_c2_nuclei'\n",
    "n_neighbors=25\n",
    "s_slide = s_sample.split('_')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#raw umap\n",
    "resolution=0.2\n",
    "results_file = f'{codedir}/Data/{s_sample}_norm_{n_neighbors}_leiden{resolution}.h5ad'  # the file that will store the analysis results\n",
    "print(results_file)\n",
    "#load leiden\n",
    "adata = sc.read_h5ad(results_file)\n",
    "rawfilename = f'{codedir}/Data/Quench/features_{s_sample}_MeanIntensity.csv'\n",
    "print(rawfilename)\n",
    "#load raw\n",
    "rdata = sc.read_csv(rawfilename)\n",
    "if s_slide == '4165NPanc':\n",
    "    rdata.obs['tissue'] = [item.split('_scene')[0].split('4165')[1] for item in rdata.obs.index]\n",
    "else:\n",
    "    rdata.obs['tissue'] = [item.split('_scene')[0] for item in rdata.obs.index]\n",
    "rdata.raw = rdata\n",
    "#log transform, reduce dimensionality (PCA)\n",
    "sc.pp.log1p(rdata)\n",
    "\n",
    "#don't remove batch effect\n",
    "#sc.pp.combat(adata,key='tissue')\n",
    "#reduce dimensionality\n",
    "sc.tl.pca(rdata, svd_solver='auto')\n",
    "'''\n",
    "sc.pp.neighbors(rdata, n_neighbors=10)\n",
    "#sc.tl.umap(rdata)\n",
    "#umap plot\n",
    "figname = f'{codedir}/Figures/{s_slide}_UmapRaw.png'\n",
    "#umap plot\n",
    "fig, ax = plt.subplots(dpi=150)\n",
    "sc.pl.umap(rdata, color='tissue',ax=ax,title='Raw')\n",
    "plt.show(fig)\n",
    "plt.tight_layout()\n",
    "fig.savefig(figname, dpi=150)\n",
    "'''\n",
    "#save\n",
    "df_raw = pd.DataFrame(rdata.raw.X, columns = adata.var.index, index = adata.obs.index)\n",
    "df_raw.merge(adata.obs,left_index=True,right_index=True).to_csv(f'{codedir}/Data/{s_sample}_raw_{n_neighbors}_leiden{resolution}.csv')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot Clustering Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load leiden clusters\n",
    "resolution = 0.2\n",
    "n_neighbors = 25\n",
    "s_sample = '808L2_c2_nuclei'\n",
    "s_slide = s_sample.split('_')[0] #'808L2'#'HER2A'#'4165NPanc'\n",
    "if s_slide == '4165NPanc':\n",
    "    adata = sc.read_h5ad(f'{codedir}/Data/{s_slide}_AF_norm_leiden.h5ad')\n",
    "else:\n",
    "    adata = sc.read_h5ad(f'{codedir}/Data/{s_sample}_norm_{n_neighbors}_leiden{resolution}.h5ad')\n",
    "adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# umap plus clustering result \n",
    "sc.set_figure_params(fontsize=16)\n",
    "figname = f'{codedir}/Figures/{s_sample}_UmapLeiden_neigh{n_neighbors}.png'\n",
    "fig, ax = plt.subplots(figsize = (4,4),dpi=300)\n",
    "sc.pl.umap(adata, color='leiden',wspace=0.2,title='Clustering w/ Leiden algorithm',ax=ax) #size=1,\n",
    "plt.tight_layout()\n",
    "plt.show(fig)\n",
    "fig.savefig(figname,bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inspect clusters\n",
    "if s_slide == '4165NPanc':\n",
    "    marker_genes = ['R1c2_nuc','R2c2_nuc','R3c2_nuc','R4c2_nuc','R5c2_nuc']\n",
    "    marker_genes_cyt = ['R1c2_cyt','R2c2_cyt','R3c2_cyt','R4c2_cyt','R5c2_cyt']\n",
    "    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,6))\n",
    "    ax1_dict = sc.pl.matrixplot(adata, marker_genes, groupby='leiden', ax=ax1,show=False,dendrogram=False, cmap='viridis', use_raw=True, colorbar_title='Nuclear \\nintensity (log)')\n",
    "    ax2_dict = sc.pl.matrixplot(adata, marker_genes_cyt, groupby='leiden',ax=ax2,show=False, dendrogram=False, cmap='viridis', use_raw=True, colorbar_title='Cytoplasm \\nintensity (log)')\n",
    "    ax1_dict['mainplot_ax'].set_xlabel('Rounds of Quenching')\n",
    "    ax1_dict['mainplot_ax'].set_xticklabels(['0','1','2','3','4'])\n",
    "    ax1_dict['mainplot_ax'].set_title('Nuclear Autofluorescence')\n",
    "    ax1_dict['mainplot_ax'].set_ylabel('Cluster ID')\n",
    "    ax2_dict['mainplot_ax'].set_xlabel('Rounds of Quenching')\n",
    "    ax2_dict['mainplot_ax'].set_xticklabels(['0','1','2','3','4'])\n",
    "    ax2_dict['mainplot_ax'].set_title('Cytoplasmic Autofluorescence')\n",
    "    fig.suptitle(f'Mean Intensity of Each Cluster, by round, location')\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{codedir}/Figures/{s_slide}_MatrixPlot_Leiden.png',dpi=200)\n",
    "else:\n",
    "    if s_slide == 'HER2A':\n",
    "        marker_genes = ['R00c2_nuclei', 'R0c2_nuclei', 'R1c2_nuclei', 'R2c2_nuclei', 'R3c2_nuclei',\n",
    "           'R4c2_nuclei', 'R5c2_nuclei','R6c2_nuclei']\n",
    "    elif s_slide == '808L2':\n",
    "        marker_genes = ['R00c2_nuclei', 'R0c2_nuclei', 'R1c2_nuclei', 'R2c2_nuclei',\n",
    "         'R3c2_nuclei', 'R4c2_nuclei', 'R5c2_nuclei', 'R6c2_nuclei',\n",
    "         'R7c2_nuclei', 'R8c2_nuclei', 'R9c2_nuclei','R10c2_nuclei']\n",
    "    fig, ax1 = plt.subplots(figsize=(6,6))\n",
    "    ax1_dict = sc.pl.matrixplot(adata, marker_genes, groupby='leiden', ax=ax1,show=False,dendrogram=False,\n",
    "                                cmap='viridis', colorbar_title='Nuclear \\nintensity (log)',\n",
    "                               vmin=5,vmax=9)\n",
    "    ax1_dict['mainplot_ax'].set_xlabel('Rounds of Quenching')\n",
    "    #ax1_dict['mainplot_ax'].set_xticklabels(['0','1','2','3','4'])\n",
    "    ax1_dict['mainplot_ax'].set_title('Nuclear Autofluorescence')\n",
    "    ax1_dict['mainplot_ax'].set_ylabel('Cluster ID')   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax =plt.subplots(figsize=(6,5))\n",
    "if s_slide == '4165NPanc':\n",
    "    ls_drop = ['8','9','10']\n",
    "    adata_new = adata[~adata.obs['leiden'].isin(ls_drop),:]\n",
    "    axes_dict = sc.pl.matrixplot(adata_new, marker_genes, groupby='leiden', ax=ax,show=False,dendrogram=False, cmap='viridis', use_raw=True, colorbar_title='Nuclear \\nintensity (log)')\n",
    "    axes_dict['mainplot_ax'].set_title('Cluster Autofluorescence')\n",
    "    axes_dict['mainplot_ax'].set_xlabel('Rounds of Quenching')\n",
    "    axes_dict['mainplot_ax'].set_ylabel('Cluster ID')\n",
    "    axes_dict['mainplot_ax'].set_xticklabels(['0','1','2','3','4'])\n",
    "    fig.savefig(f'{codedir}/Figures/{s_slide}_MatrixPlot7_Leiden.png',dpi=200)\n",
    "elif s_slide == 'HER2A':\n",
    "    ls_drop = ['6','7','8','9','10']\n",
    "    ls_scene = ['scene006','scene004','scene005','scene009']\n",
    "    adata_new = adata[((~adata.obs['leiden'].isin(ls_drop)) & (adata.obs['scene'].isin(ls_scene))),:]\n",
    "    axes_dict = sc.pl.matrixplot(adata_new, marker_genes, groupby='leiden', ax=ax,show=False,dendrogram=False, cmap='viridis', use_raw=True, colorbar_title='Nuclear \\nintensity (log)')\n",
    "    axes_dict['mainplot_ax'].set_title('Cluster Autofluorescence')\n",
    "    axes_dict['mainplot_ax'].set_xlabel('Rounds of Quenching')\n",
    "    axes_dict['mainplot_ax'].set_ylabel('Cluster ID')\n",
    "    axes_dict['mainplot_ax'].set_xticklabels(['0','1','2','3','4','5','6'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_out = adata_new.to_df()\n",
    "df_out['leiden'] = adata_new.obs.leiden\n",
    "df_out['scene'] = adata_new.obs.scene\n",
    "df_out['tissue'] = adata_new.obs.tissue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_clusters(df_cluster,df_xy):\n",
    "    s_type = df_cluster.columns[df_cluster.dtypes=='int64'][0]\n",
    "    print(s_type)\n",
    "    ls_scene = sorted(set(df_cluster.slide_scene))\n",
    "    ls_color = sorted(set(df_cluster.loc[:,s_type].dropna()))\n",
    "    ls_fig = []\n",
    "    for s_scene in ls_scene:\n",
    "        #negative cells = all cells even before dapi filtering\n",
    "        df_neg = df_xy[(df_xy.slide_scene==s_scene)]\n",
    "        #plot\n",
    "        fig, ax = plt.subplots(len(ls_color)//2 +1, 2, figsize=(12,16),dpi=200)\n",
    "        ax = ax.ravel()\n",
    "        for ax_num, s_color in enumerate(ls_color):\n",
    "            s_marker = s_color\n",
    "            #positive cells = poitive cells based on threshold\n",
    "            ls_pos_index = (df_cluster[df_cluster.loc[:,s_type]==s_color]).index\n",
    "            df_color_pos = df_neg[df_neg.index.isin(ls_pos_index)]\n",
    "            if len(df_color_pos)>=1:\n",
    "                #plot negative cells\n",
    "                ax[ax_num].scatter(data=df_neg,x='DAPI_X',y='DAPI_Y',color='silver',s=1)\n",
    "                #plot positive cells\n",
    "                ax[ax_num].scatter(data=df_color_pos, x='DAPI_X',y='DAPI_Y',color='DarkBlue',s=.5)\n",
    "                  \n",
    "                ax[ax_num].axis('equal')\n",
    "                ax[ax_num].set_ylim(ax[ax_num].get_ylim()[::-1])\n",
    "                ax[ax_num].set_xticklabels('')\n",
    "                ax[ax_num].set_yticklabels('')\n",
    "                ax[ax_num].set_title(f'{s_color} ({len(df_color_pos)} cells)')\n",
    "            else:\n",
    "                ax[ax_num].set_xticklabels('')\n",
    "                ax[ax_num].set_yticklabels('')\n",
    "                ax[ax_num].set_title(f'{s_color}  ({(0)} cells')\n",
    "        \n",
    "        fig.suptitle(s_scene)\n",
    "        ls_fig.append(fig)\n",
    "    return(ls_fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_slide = 'HER2A' #'4165NPanc'\n",
    "df_xy = pd.read_csv(f'{codedir}/Data/Quench/features_{s_slide}_CentroidXY.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load\n",
    "if s_slide == '4165NPanc':\n",
    "    df_lei = pd.read_csv(f'{codedir}/Data/{s_slide}_AF_norm_leiden.csv', index_col=0)\n",
    "else:\n",
    "    df_lei = adata.to_df()\n",
    "    df_lei['leiden'] = [int(item) for item in adata.obs.leiden]\n",
    "    df_lei['tissue'] = adata.obs.tissue\n",
    "df_lei['slide_scene'] = [item.split('_cell')[0] for item in df_lei.index]\n",
    "i_scene = len(df_lei.slide_scene.unique())/len(df_lei.tissue.unique())\n",
    "ls_fig = plot_clusters(df_lei,df_xy)\n",
    "for idx, fig in enumerate(ls_fig):\n",
    "    print(idx)\n",
    "    fig.savefig(f'{codedir}/Figures/SpatialPlots_{ls_sample[int(idx//i_scene)]}_{int(idx)}_Leiden.png')\n",
    "    plt.close(fig)\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#if not os.path.exists(f'{codedir}/Data/{s_slide}_AF_norm_leiden.csv'):\n",
    "#    df_lei.to_csv(f'{codedir}/Data/{s_slide}_AF_norm_leiden.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#histograms\n",
    "df_plot = df_cyto\n",
    "d_result = {}\n",
    "for s_channel in ['c2','c3','c4','c5']:\n",
    "    #plot single cell mean intensity\n",
    "    df_channel = df_plot.loc[:,df_plot.columns.str.contains(f'{s_channel}_')]\n",
    "    fig,ax = plt.subplots(6,1,figsize=(6,14),dpi=200)\n",
    "    ax=ax.ravel()\n",
    "    for idxs, s_sample in enumerate(ls_sample):\n",
    "        df_sample = df_channel[df_channel.index.str.contains(s_sample)]\n",
    "        df_sample_diff = pd.DataFrame(index=df_sample.index)        \n",
    "        for idx, s_round in enumerate(['R1','R2','R3','R4','R5','R6']):\n",
    "            df_round = df_sample.loc[:,df_sample.columns.str.contains(f'{s_round}c')]\n",
    "            if idx == 0:\n",
    "                i_max = np.quantile(df_round,.99)\n",
    "            sns.distplot(df_round,bins=500,kde=False,ax=ax[idx],label=s_sample)\n",
    "            i_mean = df_round.mean() \n",
    "            ax[idx].axvline(i_mean[0],ymin=0,ymax=1,color=f'C{idxs}',linestyle='--')\n",
    "            ax[idx].set_xlabel('')  \n",
    "            ax[idx].set_title(f'{s_round} {s_channel}')\n",
    "            ax[idx].set_ylabel('No. of Cells')\n",
    "            ax[idx].set_xlim(0,i_max)\n",
    "            if s_round != 'R1':\n",
    "                df_sample_diff[f'{s_round} - {s_last}'] = df_round.values - df_last.values\n",
    "            df_last = df_round\n",
    "            s_last = s_round\n",
    "        ax[idx].legend()  \n",
    "        ax[idx].set_xlabel(f'Intensity')\n",
    "        d_result.update({s_sample:df_sample_diff})\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{codedir}/Figures/4165NPanc/SingleCellDist_{s_channel}.png')\n",
    "    if s_channel != 'c2':\n",
    "        plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3.9.5",
   "language": "python",
   "name": "base"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}