20240520_NP-DCIS_cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster TMAs\n",
    "\n",
    "**Samples:** \n",
    "- NP8\n",
    "- NP9\n",
    "\n",
    "**Method**: \n",
    "\n",
    "- threshold and gate cell types\n",
    "\n",
    "- marker expression in tumor/stroma\n",
    "\n",
    "- leiden cluster on good biomarkers for celltypes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load libraries\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import seaborn as sns\n",
    "import scipy\n",
    "import scanpy as sc\n",
    "import matplotlib as mpl\n",
    "from matplotlib import cm\n",
    "import math\n",
    "mpl.rc('figure', max_open_warning = 0)\n",
    "codedir = '/home/groups/graylab_share/Chin_Lab/ChinData/engje/Data/20200000/20200406_JP-TMAs'#os.getcwd()\n",
    "import phenograph\n",
    "os.chdir('/home/groups/graylab_share/Chin_Lab/ChinData/engje/Data')\n",
    "from mplex_image import visualize as viz, process, preprocess, normalize, mics, mpimage\n",
    "np.random.seed(126)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "codedir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change to correct directory\n",
    "datadir = f'/home/groups/graylab_share/Chin_Lab/ChinData/engje/Data/20200000/20200406_JP-TMAs/NP_DCIS'\n",
    "\n",
    "s_date = '20240520'\n",
    "os.chdir(codedir)\n",
    "if not os.path.exists(s_date):\n",
    "    os.mkdir(s_date)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of contents <a name=\"contents\"></a>\n",
    "2. [load data](#cluster)\n",
    "3. [select rois](#load)\n",
    "5. [single cell Umap](#l7) **start here**\n",
    "6. [Leiden celltypes](#l8)\n",
    "[annotate leiden celltypes](#lbar)\n",
    "7. [Tissue means](#tissue)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NP Tissues  <a name=\"cluster\"></a>\n",
    "\n",
    "\n",
    "filter markers, select ROIs\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.listdir(datadir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mean int\n",
    "os.chdir(datadir)\n",
    "ls_sample = ['NP9-DCIS2','NP8-DCIS','NP8-IDC','NP9-IDC2'] \n",
    "df_mi = pd.DataFrame()\n",
    "df_xy = pd.DataFrame()\n",
    "for s_sample in ls_sample:\n",
    "    df_mi2 = pd.read_csv(f'{datadir}/features_{s_sample}_FilteredMeanIntensity.txt',\n",
    "                         index_col=0, sep='\\t')\n",
    "    #df_mi2['slide_scene'] = [item.split('_cell')[0] for item in df_mi2.index]\n",
    "    df_xy2 = pd.read_csv(f'{datadir}/features_{s_sample}_CentroidXY.csv',index_col=0)\n",
    "    df_xy2['slide_scene'] = [item.split('_cell')[0] for item in df_xy2.index]\n",
    "    df_mi =pd.concat([df_mi,df_mi2])\n",
    "    df_xy =pd.concat([df_xy,df_xy2])\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df_mi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#shorten and check\n",
    "df_mi.columns = [item.split('_')[0] for item in df_mi.columns]\n",
    "print(len(df_mi.columns))\n",
    "df_mi.columns[df_mi.columns.duplicated()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi = df_mi.merge(df_xy,left_index=True, right_index=True)\n",
    "#df_mi.rename({'nuclei_area':'area','nuclei_eccentricity':'eccentricity'},axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save\n",
    "#s_out = f'{codedir}/data/20211221_JP-TMA2_FilteredMeanIntensity.csv'\n",
    "s_out = f'{codedir}/data/20240520_NP-Tissues_FilteredMeanIntensity.csv'\n",
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_mi.to_csv(s_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ROI selection <a name=\"load\"></a>\n",
    "\n",
    "select virtual cores 1.2mm diameter\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_out = f'{codedir}/data/20240520_NP-Tissues_FilteredMeanIntensity.csv'\n",
    "df_mi = pd.read_csv(s_out,index_col=0)\n",
    "df_mi['slide'] = [item.split('_')[0] for item in df_mi.index]\n",
    "df_mi['CK19_Pos'] = df_mi.CK19>1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## visualize\n",
    "# a 1.2 um core\n",
    "radius = 1200*3.25/2 #3.25  pixel per um\n",
    "\n",
    "def plot_pos(df_pos,ls_color,df_xy):\n",
    "    #plot\n",
    "    fig, ax = plt.subplots(figsize=(5,5)) #figsize=(18,12)\n",
    "    #plot negative cells\n",
    "    ax.scatter(data=df_xy,x='DAPI_X',y='DAPI_Y',color='silver',s=1)\n",
    "    for idx,s_color in enumerate(ls_color):\n",
    "        #positive cells = positive cells based on threshold\n",
    "        ls_pos_index = (df_pos[df_pos.loc[:,s_color]]).index\n",
    "        df_color_pos = df_xy[df_xy.index.isin(ls_pos_index)]\n",
    "        if len(df_color_pos)>=1:\n",
    "            #plot positive cells\n",
    "            ax.scatter(data=df_color_pos, x='DAPI_X',y='DAPI_Y',color=f'C{idx}',s=.25,alpha=0.8)\n",
    "    ax.axis('equal')\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    ax.set_title(f'{s_color}')\n",
    "    fig.suptitle(s_scene)\n",
    "    return(fig, ax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load hand drawn ROIS\n",
    "df_roi = pd.DataFrame()\n",
    "for s_file in os.listdir(f'{codedir}/NP9-ROI_selection'):\n",
    "    if s_file.find('exclude') >-1:\n",
    "        df = pd.read_csv(f'{codedir}/NP9-ROI_selection/{s_file}',index_col=0)\n",
    "        s_roi = s_file.split('_')[-1].split('.')[0]\n",
    "        df['ROI'] = s_roi\n",
    "        df_roi = pd.concat([df_roi,df])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted(df_mi.slide_scene.unique())\n",
    "d_type = {#'NP8-DCIS_scene001':'skip', #DCIS\n",
    " #'NP8-DCIS_scene002':'DCIS',\n",
    " #'NP8-DCIS_scene003':'DCIS',\n",
    " #'NP8-DCIS_scene004':'DCIS',\n",
    " #'NP8-DCIS_scene005':'DCIS',\n",
    " #'NP8-DCIS_scene006':'skip', #DCIS\n",
    " #'NP8-DCIS_scene007':'skip', #DCIS\n",
    " #'NP8-DCIS_scene008':'DCIS',\n",
    " #'NP8-DCIS_scene009':'DCIS',\n",
    " #'NP8-DCIS_scene010':'DCIS',\n",
    " #'NP8-IDC_scene001':'IDC',\n",
    " #'NP8-IDC_scene002':'IDC',\n",
    "# 'NP8-IDC_scene003':'IDC',\n",
    "# 'NP8-IDC_scene004':'skip', #IDC\n",
    "# 'NP8-IDC_scene005':'skip', #IDC\n",
    " #'NP9-DCIS2_scene001':'Normal',\n",
    " #'NP9-DCIS2_scene002':'DCIS',\n",
    " #'NP9-DCIS2_scene003':'DCIS',\n",
    " 'NP9-IDC2_scene001':'Normal',\n",
    " #'NP9-IDC2_scene002':'DCIS',\n",
    " 'NP9-IDC2_scene003':'IDC'}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot CK19\n",
    "for s_scene, s_type in d_type.items():\n",
    "    df_scene = df_mi[df_mi.slide_scene==s_scene]\n",
    "    fig, ax = plot_pos(df_mi,['CK19_Pos'],df_scene)\n",
    "    ax.set_title(f'{s_type} CK19+')\n",
    "    break\n",
    "# hand drawn rois for that scene\n",
    "df_roi_scene = df_roi[df_roi.index.str.contains(s_scene)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_scene = df_mi[df_mi.slide_scene==s_scene].copy()\n",
    "for s_roi in df_roi_scene.ROI.unique():\n",
    "    ls_index = df_roi_scene[df_roi_scene.ROI==s_roi].index\n",
    "    df_scene[s_roi] = df_scene.index.isin(ls_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #select circular ROIs\n",
    "# dd_center= { #'NP8-DCIS_scene001':{'DCIS1':(10000,7000)},\n",
    "#     'NP8-IDC_scene001':{'Tumor1':(14000,8000),'Border1':(11000,3500),'Border2':(7000,8000)},\n",
    "#      'NP8-IDC_scene002':{'Tumor1':(13500,7000),'Tumor2':(17000,4000),'Border1':(9500,8000)}, \n",
    "#     'NP8-IDC_scene003':{'Tumor1':(6000,5000),'Tumor2':(12500,7500),'Border1':(2000,6500),'Border2':(17000,6000)},\n",
    "#     'NP9-DCIS2_scene001':{'Normal1':(2000,2000),'Normal2':(2500,7000),'Normal3':(11000,4000),'Normal4':(16000,5000)},\n",
    "#     'NP9-IDC2_scene001':{'Normal1':(5500,8500),'Normal2':(3000,11700)},\n",
    "#     'NP9-IDC2_scene003':{'Tumor1':(9500,18500),'Border1':(10500,14500)},\n",
    "# }\n",
    "# #slow #select circular ROIs\n",
    "# dd_result = {}\n",
    "# for s_scene, d_center in dd_center.items():\n",
    "#     d_result={}\n",
    "#     df_scene = df_mi[df_mi.slide_scene==s_scene].copy()\n",
    "#     for idx, (s_roi, center) in enumerate(d_center.items()):\n",
    "#         # Calculate Euclidean distance\n",
    "#         for (index, Series) in df_scene.iterrows():\n",
    "#             p = [Series.DAPI_X,Series.DAPI_Y]\n",
    "#             dist =  math.dist(p, center)\n",
    "#             df_scene.loc[index,'distance'] = dist\n",
    "#         df_scene[f'ROI{idx}'] = df_scene.distance < radius\n",
    "#         ls_index = df_scene[df_scene.loc[:,f'ROI{idx}']].index.tolist()\n",
    "#         d_result.update({s_roi:ls_index})\n",
    "#     dd_result.update({s_scene:d_result})\n",
    "#     #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot ROIS\n",
    "fig, ax = plot_pos(df_scene,['CK19_Pos','Normal1','Normal2','Normal3','Normal4'],df_scene) #,'ROI2','ROI3','ROI0','ROI1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #save json\n",
    "# import json\n",
    "# with open('NP_ROI_data.json', 'w') as f:\n",
    "#     json.dump(dd_result, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #load json\n",
    "# import json\n",
    "\n",
    "# with open('NP_ROI_data.json') as f:\n",
    "#     dd_result = json.load(f)\n",
    "#     print(dd_result.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #filter\n",
    "# df_mi['ROI'] = pd.NA\n",
    "# for s_scene, d_result in dd_result.items():\n",
    "#     for s_roi, ls_index in d_result.items():\n",
    "#         df_mi.loc[ls_index,'ROI'] = f'{s_scene}_{s_roi}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## filter hand drawn\n",
    "df_mi['ROI'] = pd.NA\n",
    "for s_roi in df_roi.ROI.unique():\n",
    "    print(s_roi)\n",
    "    ls_index = df_roi[df_roi.ROI==s_roi].index\n",
    "    s_scene = ls_index[0].split('_cell')[0]\n",
    "    df_mi.loc[ls_index,'ROI'] = f'{s_scene}_{s_roi}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save ROIS\n",
    "s_out = f'{codedir}/data/20240521_NP-Tissues-ROIs_FilteredMeanIntensity.csv'\n",
    "s_out = f'{codedir}/data/20240522_NP9-IDC2-ROIs_FilteredMeanIntensity.csv'\n",
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_mi[df_mi.ROI.notna()].to_csv(s_out)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Umap <a name=\"l7\"></a>\n",
    "\n",
    "umap projection of single cells\n",
    "\n",
    "first use all markers, all subtypes, just to evaluate markers\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scanpy as sc\n",
    "import os\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(datadir)\n",
    "#load ROIS\n",
    "s_out = f'{codedir}/data/20240521_NP-Tissues-ROIs_FilteredMeanIntensity.csv'\n",
    "s_out = f'{codedir}/data/20240522_NP9-IDC2-ROIs_FilteredMeanIntensity.csv'\n",
    "df_norm  = pd.read_csv(f'{s_out}',index_col=0) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ER+\n",
    "df_exclude = pd.DataFrame() #pd.read_csv('exclude_JP-TMA-1_ER+_34markers_leiden0.6_20.csv',index_col=0)\n",
    "\n",
    "ls_exclude = []\n",
    "#ls_exclude_all = sorted(set(df_norm[df_norm.loc[:,'PD1'] > 1365].loc[:,'PD1'].index.tolist()).union(set(ls_exclude))) #PD1\n",
    "#df_exclude = df_exclude.append(pd.DataFrame(index=ls_exclude_all,columns=df_exclude.columns))\n",
    "df_exclude_er = df_exclude.copy()\n",
    "print(len(df_exclude_er))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#32 features\n",
    "#[ 'NP8-IDC' 'NP9-IDC2']\n",
    "s_sample = s_out.split('data/')[1].split('_Filtered')[0]\n",
    "s_type = 'NP9-IDC2'#'NP8-IDC'# 'all'\n",
    "d_roi = dict(zip(df_mi.index,df_mi.ROI))\n",
    "#almost same features as IMC\n",
    "ls_col = [ 'CD31', 'CD44', 'CD45','CD20','CD4','CD8',#'CD3',\n",
    "        'CK14','CK5','EGFR', 'CK19',  'CK8','CK7','Ecad',  #\n",
    "       'Ki67','CD68', 'aSMA','ColI','HER2', 'ER','PgR', #np8 is TNBC?\n",
    "        'Vim'] #,'area'\n",
    "print(len(ls_col))\n",
    "adata = sc.AnnData(df_norm.loc[(~df_norm.index.isin(df_exclude.index)),ls_col]) \n",
    "adata.obs['slide_scene'] = [item.split('_cell')[0] for item in adata.obs.index]\n",
    "adata.obs['ROI'] = adata.obs.index.map(d_roi)\n",
    "adata.obs['subtype'] =[item.split('_')[-1][0:-1] for item in adata.obs.ROI]\n",
    "adata.obs['TMA'] = [item.split('_')[0] for item in adata.obs.index]\n",
    "# #two subtypes\n",
    "if s_type == 'all':\n",
    "    adata = adata[adata.obs.subtype.isin(['Normal', 'Border', 'Tumor'])]\n",
    "else:\n",
    "    adata = adata[adata.obs.TMA.isin([s_type])]  \n",
    "print(len(adata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include two TNBC TMAs\n",
    "print(adata.obs['TMA'].unique())\n",
    "print(adata.obs['subtype'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#raw\n",
    "b_scale = True\n",
    "adata.raw = adata\n",
    "#reduce dimensionality\n",
    "sc.tl.pca(adata, svd_solver='auto')\n",
    "fig,ax=plt.subplots(figsize=(3.5,5))\n",
    "sc.pl.highest_expr_genes(adata, n_top=48,ax=ax,save=f'Expression_{len(ls_col)}.png')\n",
    "plt.tight_layout()\n",
    "sc.pl.pca_variance_ratio(adata,n_pcs=32, log=True)\n",
    "\n",
    "if b_scale:\n",
    "    sc.pp.scale(adata, zero_center=False, max_value=20)\n",
    "    s_sample = s_sample + '_s'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(ls_col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#umap\n",
    "for n_neighbors in [30,15]:\n",
    "    results_file = f'{s_sample}_{n_neighbors}neighbors_{s_type}_{len(ls_col)}markers.h5ad'\n",
    "    if not os.path.exists(results_file):\n",
    "        print('calc umap')\n",
    "        # calculate neighbors \n",
    "        sc.pp.neighbors(adata, n_neighbors=n_neighbors) #, method='rapids'\n",
    "        sc.tl.umap(adata)\n",
    "    else:\n",
    "        print('loading umap')\n",
    "        adata = sc.read_h5ad(results_file)\n",
    "        adata.obs['slide_scene'] = [item.split('_cell')[0] for item in adata.obs.index]\n",
    "        adata.obs['subtype'] = adata.obs.slide_scene.map(d_cyc_sub).fillna('other').replace(d_replace)\n",
    "        adata.obs['celltype'] = df_man.loc[adata.obs.index,'celltype']\n",
    "    #save results\n",
    "\n",
    "    #color by markers\n",
    "    #figname = f\"Umap{s_sample}_markers_{n_neighbors}neighbors_{len(ls_col)}markers.png\"\n",
    "    #sc.pl.umap(adata, color=ls_col,vmin='p1.5',vmax='p98.5',save=figname,ncols=6)\n",
    "    #color by TMA\n",
    "    figname = f\"Umap{s_sample}_TMA_{n_neighbors}neighbors_{s_type}_{len(ls_col)}markers.png\"\n",
    "    fig,ax = plt.subplots(figsize=(3,2), dpi=200)\n",
    "    sc.pl.umap(adata, color='TMA',wspace=.25,save=figname,ax=ax)\n",
    "    \n",
    "    figname = f\"Umap{s_sample}_subtype_{n_neighbors}neighbors_{s_type}_{len(ls_col)}markers.png\"\n",
    "    fig,ax = plt.subplots(figsize=(3,2), dpi=200)\n",
    "    sc.pl.umap(adata, color='subtype',wspace=.25,save=figname,ax=ax)\n",
    "    \n",
    "    if not os.path.exists(results_file):\n",
    "        print(\"writing\")\n",
    "        adata.write(results_file)\n",
    "    break\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Cluster Leiden <a name=\"l8\"></a>\n",
    "\n",
    "cluster umap projection of single cells\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#cluster\n",
    "#s_sample = '20240521_NP-Tissues-ROIs_s'#'\n",
    "#n_neighbors = 30 #\n",
    "n_markers = len(ls_col) #19 #\n",
    "#s_type = 'NP9-IDC2'#'NP8-IDC'#'all' #\n",
    "resolution = 0.6\n",
    "results_file = f'{s_sample}_{n_neighbors}neighbors_{n_markers}markers_{s_type}_leiden{resolution}.h5ad'\n",
    "if not os.path.exists(results_file):\n",
    "    print('clustering')\n",
    "    sc.tl.leiden(adata,resolution=resolution)\n",
    "    adata.write(results_file)\n",
    "else:\n",
    "    print('loading leiden')\n",
    "    adata = sc.read_h5ad(results_file)\n",
    "\n",
    "\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}_nolegend.png'\n",
    "sc.pl.umap(adata, color='leiden',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '),legend_loc=None)\n",
    "\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "sc.pl.umap(adata, color='leiden',ax=ax)\n",
    "\n",
    "figname = f\"Umap{s_sample}_subtype_{n_neighbors}neighbors_{s_type}_{len(ls_col)}markers.png\"\n",
    "fig,ax = plt.subplots(figsize=(3,2), dpi=200)\n",
    "sc.pl.umap(adata, color='subtype',wspace=.25,save=figname,ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# sc.pl.umap(adata, color='Vim',vmin='p1.5',vmax='p98.5')\n",
    "# sc.pl.umap(adata, color='ColI',vmin='p1.5',vmax='p98.5')\n",
    "# sc.pl.umap(adata, color='Ecad',vmin='p1.5',vmax='p98.5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#legend\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "sc.pl.umap(adata, color='leiden',ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annotate Leiden <a name=\"lbar\"></a>\n",
    "\n",
    "annotate epitelial, immune a stroma cell types\n",
    "\n",
    "also visualize cluster results on tissue scatter plots. Additional visualization done using 20201018_JP-TMAs_napari.py,\n",
    "to visualize ome.tiff created in 20201005_JP-TMA_Pipeline.py\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(f'{datadir}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load ROIS\n",
    "s_out = f'{codedir}/data/20240521_NP-Tissues-ROIs_FilteredMeanIntensity.csv'\n",
    "#s_out = f'{codedir}/data/20240522_NP9-IDC2-ROIs_FilteredMeanIntensity.csv'\n",
    "df_norm  = pd.read_csv(f'{s_out}',index_col=0) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load leiden\n",
    "s_sample = s_out.split('_Filtered')[0].split('data/')[1] + '_s' #'20240521_NP-Tissues-ROIs_s'#''20240522_NP9-IDC2-ROIs_s'#\n",
    "n_neighbors = 30 #\n",
    "resolution = 0.6\n",
    "if s_sample == '20240521_NP-Tissues-ROIs_s':\n",
    "    s_type = 'NP8-IDC'#\n",
    "    n_markers =19# 21 #\n",
    "else:\n",
    "    s_type='NP9-IDC2'#'NP9-IDC2'#'all' #\n",
    "    n_markers = 21 #\n",
    "results_file = f'{s_sample}_{n_neighbors}neighbors_{n_markers}markers_{s_type}_leiden{resolution}.h5ad'\n",
    "adata = sc.read_h5ad(results_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# #unmarked stroma\n",
    "# #df_str = pd.read_csv(f'annotated_JP-TMA1-1_unmarked_stroma_{s_type}.csv',index_col=0)\n",
    "# #df_str = pd.read_csv(f'annotated_JP-TMAs_unmarked_stroma_{s_type}.csv',index_col=0)\n",
    "# fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "# figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}_legend.png'\n",
    "# sc.pl.umap(adata, color='leiden',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '))\n",
    "\n",
    "# results_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ls_col = ['CD31','Ecad','ColI','Vim','CD45']\n",
    "# figname = f\"Umap_markers_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.png\"\n",
    "# title=figname.split('.png')[0].replace('_',' ')\n",
    "# sc.pl.umap(adata, color=ls_col,vmin='p1.5',vmax='p99',ncols=5,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# annotate\n",
    "df_p = adata.to_df()\n",
    "df_p['leiden'] = adata.obs['leiden'].astype('str')\n",
    "\n",
    "if results_file =='20240521_NP-Tissues-ROIs_s_30neighbors_21markers_NP9-IDC2_leiden0.6.h5ad':\n",
    "    ls_drop = ['13','11','12'] #13 missing rounds, 3, 7 is AF\n",
    "    df_p.loc[(df_p.leiden=='4') & (df_p.CD20>9),'leiden'] = '4b'\n",
    "    df_p.loc[(df_p.leiden=='4') & (df_p.CD4>5.5),'leiden'] = '4t'\n",
    "    tum_clust = ['5','0','2']\n",
    "    str_clust = ['8','3','11','12','7']\n",
    "    imm_clust = ['4','4t','4b'] \n",
    "    endo_clust = ['6']\n",
    "    fb_clust = ['9','10','1']\n",
    "    d_named = {'3':'AF','7':'AF','5':'Normal Breast','0':'ER+ tum.','2':'ER+ tum.','6':'endothelial','4':'macrophage',\n",
    "              '4t':'CD4 T cell','4b':'B cell','9':'Vim+ FB','10':'ColI+ FB','1':'FB','8':'Quies. Str.',\n",
    "              }\n",
    "elif results_file =='20240521_NP-Tissues-ROIs_s_30neighbors_21markers_NP8-IDC_leiden0.6.h5ad':\n",
    "    ls_drop = ['16','14','7','13','12']#11 take out PGR/ER #16 missing rounds, 14 and 13 stiching error, 7 AF including necrotic, 12 necrotic\n",
    "    tum_clust = ['15']\n",
    "    str_clust = [] \n",
    "    imm_clust = [] \n",
    "    endo_clust = []\n",
    "    fb_clust = []\n",
    "    d_named = {'15':'myoep','7':'AF','12':'necrotic'}\n",
    "elif results_file =='20240521_NP-Tissues-ROIs_s_30neighbors_19markers_NP8-IDC_leiden0.6.h5ad':\n",
    "    print('used')\n",
    "    ls_drop = ['14','12','11','8'] #14 missing rounds, '12' unfortunate CD20 BG 11  AF including necrotic, 8 necrotic\n",
    "    #drop scene003 tumor 1 clust 4- weird bg\n",
    "    df_p.loc[(df_p.leiden=='6') & (df_p.CD20>11.5),'leiden'] = '6b'\n",
    "    df_p.loc[(df_p.leiden=='6') & (df_p.CD4>7),'leiden'] = '6t'\n",
    "    #df_p.loc[(df_p.leiden=='10') & (df_p.CD31>4),'leiden'] = '10e'\n",
    "    tum_clust = ['13','9','3','0','1','4']\n",
    "    str_clust = ['2'] \n",
    "    imm_clust = ['7','6','6b','6t'] \n",
    "    endo_clust = ['10']\n",
    "    fb_clust = ['5']\n",
    "    d_named = {'13':'myoep','11':'AF','10':'endothelial','9':'Prolif. tum.','8':'necrotic','7':'macrophage',\n",
    "              '6':'CD8 T cell','6b':'B cell','6t':'CD4 T cell','5':'FB','3':'Luminal tum.','0':'Luminal tum.',\n",
    "              '1':'Luminal tum.','2':'Quies. Str.','4':'Luminal tum.',}\n",
    "elif results_file == '20240522_NP9-IDC2-ROIs_s_30neighbors_21markers_NP9-IDC2_leiden0.6.h5ad':\n",
    "    print('also used')\n",
    "    df_p.loc[(df_p.leiden=='3') & (df_p.CD20>9),'leiden'] = '3b'\n",
    "    ls_drop = ['7','8','15','17','16','13','14'] #7 and 8 is AF, 15 is missing rounds, 16 and 17 is AF, 14 is af, 13 floating\n",
    "    tum_clust = ['0','5','10']\n",
    "    str_clust = ['9'] \n",
    "    imm_clust = ['12','3','3b'] #12 mac\n",
    "    endo_clust = ['6']\n",
    "    fb_clust = ['1','4','2','11',]\n",
    "    d_named = {'10':'myoep','5':'Normal breast luminal','0':'ER+ Tumor','6':'endothelial','3b':'B cell','3':'CD4 T cell',\n",
    "              '12':'macrophage','11':'Vim+ FB','4':'ColI++ FB','2':'ColI+ FB','1':'FB','9':'Quies. Str.'}\n",
    "else:\n",
    "    print('new one')\n",
    "    tum_clust = []\n",
    "    str_clust = [] \n",
    "    imm_clust = [] \n",
    "    endo_clust = []\n",
    "    fb_clust = []\n",
    "\n",
    "    \n",
    "#df_p.groupby('leiden').mean().Ecad.sort_values(ascending=False).index\n",
    "#set(df_p.groupby('leiden').mean().CK19.sort_values(ascending=False).index) - set(tum_clust + imm_clust )\n",
    "#set(tum_clust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clustermap\n",
    "if s_type=='NP9-IDC2':\n",
    "    ls_drop_col = ['CD8','HER2','Ki67']\n",
    "else:\n",
    "    ls_drop_col = []\n",
    "b_annot = False\n",
    "d_replace = {}\n",
    "d_replace.update(dict(zip(tum_clust,[f'epithelial' for item in tum_clust])))\n",
    "d_replace.update(dict(zip(str_clust,[f'stromal' for item in str_clust])))\n",
    "d_replace.update(dict(zip(imm_clust,[f'immune' for item in imm_clust])))\n",
    "d_replace.update(dict(zip(endo_clust,[f'endothelial' for item in endo_clust])))\n",
    "d_replace.update(dict(zip(fb_clust,[f'fibroblast' for item in fb_clust])))\n",
    "df_plot = df_p.loc[~df_p.leiden.isin(ls_drop),~df_p.columns.isin(ls_drop_col)].groupby('leiden').mean()\n",
    "\n",
    "df_plot.index.name = f'leiden {resolution}'\n",
    "d_color = dict(zip(['endothelial', 'epithelial', 'fibroblast', 'immune', 'stromal'],sns.color_palette()[0:5]))\n",
    "if b_annot:\n",
    "    row_colors = df_plot.index.astype('str').map(d_replace).map(d_color)\n",
    "    g = sns.clustermap(df_plot,z_score=1,figsize=(9,7),cmap='viridis',row_colors=row_colors,#method='single',\n",
    "                   vmin=-1.5,vmax=1.5)\n",
    "    for label in d_color.keys():\n",
    "        g.ax_row_dendrogram.bar(0, 0, color=d_color[label],\n",
    "                                label=label, linewidth=0)\n",
    "    g.ax_row_dendrogram.legend(loc=\"right\", ncol=1)\n",
    "    g.ax_heatmap.set_ylabel('')\n",
    "    marker_genes = df_plot.iloc[:,g.dendrogram_col.reordered_ind].columns.tolist()\n",
    "    categories_order = df_plot.iloc[g.dendrogram_row.reordered_ind,:].index.tolist()\n",
    "    labels = [f\"{item}: {d_named[item]}\" for item in categories_order]\n",
    "    g.ax_heatmap.set_yticklabels(labels,rotation=0)\n",
    "else:\n",
    "    g = sns.clustermap(df_plot,z_score=1,figsize=(9,7),cmap='viridis',#row_colors=row_colors,#method='single',\n",
    "                   vmin=-1.5,vmax=1.5)\n",
    "    marker_genes = df_plot.iloc[:,g.dendrogram_col.reordered_ind].columns.tolist()\n",
    "    categories_order = df_plot.iloc[g.dendrogram_row.reordered_ind,:].index.tolist()\n",
    "g.savefig(f'{codedir}/{s_date}/{s_sample}_clustermap_{n_neighbors}_{n_markers}markers_{s_type}_leiden{resolution}.png',dpi=200)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(1.4,3.7),dpi=200)\n",
    "df_p.groupby('leiden').count().loc[categories_order[::-1]].iloc[:,1].plot(kind='barh',title='Cell No.',ax=ax,width=0.7)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/{s_date}/barplot_{s_sample}_{n_neighbors}neighbors_{len(marker_genes)}markers_leiden{resolution}_{s_type}.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "s_clust = '8'\n",
    "df_p['slide_scene'] = [item.split('_cell')[0] for item in df_p.index]\n",
    "df_p.groupby(['leiden','slide_scene']).count().loc[s_clust].iloc[:,1].sort_values(ascending=False)[0:30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mplex_image import analyze\n",
    "df_data = df_norm\n",
    "colors = mpl.cm.tab10.colors + mpl.cm.Accent.colors + mpl.cm.tab20b.colors + mpl.cm.tab20c.colors\n",
    "\n",
    "df_pos = analyze.celltype_to_bool(df_p,'leiden')\n",
    "df_pos.columns = [str(item) for item in df_pos.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#plot all groups spatially - leiden \n",
    "ls_clust = df_pos.columns.tolist()#['8','11','12','7','5','10','1','9']#'0','2',\n",
    "for s_slide in sorted(set(df_p.slide_scene)):\n",
    "    #s_slide = 'JP-TMA1-1_scene029'\n",
    "    fig,ax = plt.subplots(figsize=(5,4.5),dpi=200)\n",
    "    #plot negative cells\n",
    "    df_scene = df_data[df_data.index.str.contains(s_slide)]\n",
    "    ax.scatter(data=df_scene,x='DAPI_X',y='DAPI_Y',color='silver',s=0.1,label=f'')\n",
    "    #for idxs, s_color_int in enumerate(range(len(df_pos.columns))):\n",
    "    for idxs, s_color in enumerate(ls_clust): #enumerate(df_pos.columns):\n",
    "        #s_color = s_clust\n",
    "        s_color = str(s_color)\n",
    "        if len(df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ls_index = df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])].index\n",
    "            ax.scatter(data=df_data.loc[ls_index],x='DAPI_X',y='DAPI_Y',label=f'{s_color}',s=0.1,color=colors[idxs])\n",
    "        #break\n",
    "    ax.set_title(f\"{s_slide}\", fontsize=16) # \\n {d_a[s_slide]}\n",
    "    ax.axis('equal')\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    #ax.set_xticklabels('')\n",
    "    #ax.set_yticklabels('')\n",
    "    #break\n",
    "    plt.legend(markerscale=10,framealpha=.5,bbox_to_anchor=(1,1)) \n",
    "    #fig.savefig(f'{codedir}/{s_date}/{s_slide}_leiden{resolution}_scatterplot.png')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ls_merge = tum_clust# fb_clust #tum_clust #['0','12'] #16\n",
    "adata.obs.leiden = df_p.leiden\n",
    "adata.obs['test'] = (adata.obs.leiden.isin(ls_merge)).replace({True:1,False:0})\n",
    "sc.pl.umap(adata, color='test',title='tumor')\n",
    "print(sum(adata.obs['test'])/len(adata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adata.obs['leiden'] = df_p['leiden']\n",
    "fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}.png'\n",
    "sc.pl.umap(adata, color='leiden',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for s_merge in tum_clust:\n",
    "#     ls_merge = [s_merge]\n",
    "#     adata.obs['test'] = (adata.obs.leiden.isin(ls_merge)).replace({True:1,False:0})\n",
    "#     sc.pl.umap(adata, color='test',title=s_merge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## add celltypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_p.loc[df_p.leiden.isin(tum_clust),'leidencelltype3'] = 'epithelial'\n",
    "df_p.loc[df_p.leiden.isin(str_clust),'leidencelltype3'] = 'stromal'\n",
    "df_p.loc[df_p.leiden.isin(fb_clust),'leidencelltype3'] = 'stromal'\n",
    "df_p.loc[df_p.leiden.isin(endo_clust),'leidencelltype3'] = 'stromal'\n",
    "df_p.loc[df_p.leiden.isin(imm_clust),'leidencelltype3'] = 'immune'\n",
    "adata.obs['leidencelltype3'] = df_p.leidencelltype3\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}leidencelltype3.png'\n",
    "sc.pl.umap(adata, color='leidencelltype3',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '),legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_p.loc[df_p.leiden.isin(tum_clust),'leidencelltype5'] = 'epithelial'\n",
    "df_p.loc[df_p.leiden.isin(str_clust),'leidencelltype5'] = 'stromal'\n",
    "df_p.loc[df_p.leiden.isin(fb_clust),'leidencelltype5'] = 'fibroblast'\n",
    "df_p.loc[df_p.leiden.isin(endo_clust),'leidencelltype5'] = 'endothelial'\n",
    "df_p.loc[df_p.leiden.isin(imm_clust),'leidencelltype5'] = 'immune'\n",
    "adata.obs['leidencelltype5'] = df_p.leidencelltype5\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}_leidencelltype5.png'\n",
    "sc.pl.umap(adata, color='leidencelltype5',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#visulaize\n",
    "ls_col = ['CD31','Ecad','Vim','ColI','CD45']\n",
    "figname = f\"Umap_markers_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.png\"\n",
    "title=figname.split('.png')[0].replace('_',' ')\n",
    "sc.pl.umap(adata, color=ls_col,vmin='p1.5',vmax='p99',ncols=5,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_p.groupby('slide').leidencelltype5.value_counts(normalize=True).unstack().sort_values(by='epithelial')[25:40]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_p.loc[df_p.leiden.isin(tum_clust),'leidencelltype4'] = 'epithelial'\n",
    "df_p.loc[df_p.leiden.isin(str_clust),'leidencelltype4'] = 'stromal'\n",
    "df_p.loc[df_p.leiden.isin(fb_clust),'leidencelltype4'] = 'stromal'\n",
    "df_p.loc[df_p.leiden.isin(endo_clust),'leidencelltype4'] = 'endothelial'\n",
    "df_p.loc[df_p.leiden.isin(imm_clust),'leidencelltype4'] = 'immune'\n",
    "adata.obs['leidencelltype4'] = df_p.leidencelltype4\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}_leidencelltype4.png'\n",
    "sc.pl.umap(adata, color='leidencelltype4',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '),legend_loc='on data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.leiden.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_roi = dict(zip(adata.obs.index,adata.obs.ROI))\n",
    "df_p['ROI'] = df_p.index.map(d_roi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #drop scene003 tumor 1 clust 4- weird bg - include as tumor\n",
    "# if s_type == 'NP8-IDC':\n",
    "#     df_p = df_p[df_p.ROI!='NP8-IDC_scene003_Tumor1']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_p['annotated'] = df_p.leiden.map(d_named).fillna('AF')\n",
    "#plot all groups spatially - annotated \n",
    "colors = mpl.cm.tab10.colors + mpl.cm.Accent.colors #mpl.cm.tab20b.colors + mpl.cm.tab20c.colors\n",
    "\n",
    "df_pos = analyze.celltype_to_bool(df_p,'annotated')\n",
    "ls_clust = df_pos.columns.tolist()#['8','11','12','7','5','10','1','9']#'0','2',\n",
    "for s_slide in sorted(set(df_p.slide_scene)):\n",
    "    #s_slide = 'JP-TMA1-1_scene029'\n",
    "    fig,ax = plt.subplots(figsize=(6,6),dpi=300)\n",
    "    #plot negative cells\n",
    "    df_scene = df_data[df_data.index.str.contains(s_slide)]\n",
    "    ax.scatter(data=df_scene,x='DAPI_X',y='DAPI_Y',color='silver',s=0.1,label=f'')\n",
    "    #for idxs, s_color_int in enumerate(range(len(df_pos.columns))):\n",
    "    for idxs, s_color in enumerate(ls_clust): #enumerate(df_pos.columns):\n",
    "        #s_color = s_clust\n",
    "        s_color = str(s_color)\n",
    "        if len(df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ls_index = df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])].index\n",
    "            ax.scatter(data=df_data.loc[ls_index],x='DAPI_X',y='DAPI_Y',label=f'{s_color}',s=0.3,color=colors[idxs])\n",
    "        #break\n",
    "    ax.set_title(f\"{s_slide}\", fontsize=16) # \\n {d_a[s_slide]}\n",
    "    ax.axis('equal')\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    #ax.set_xticklabels('')\n",
    "    #ax.set_yticklabels('')\n",
    "    #break\n",
    "    plt.legend(markerscale=10,framealpha=.5,bbox_to_anchor=(1,1)) \n",
    "    #fig.savefig(f'{codedir}/{s_date}/{s_slide}_leiden{resolution}_scatterplot.png')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save\n",
    "s_out = f'{s_sample}_LeidenClustering_neighbors{n_neighbors}_resolution{resolution}_markers{n_markers}_{s_type}.csv'\n",
    "df_out = df_p\n",
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_out.to_csv(s_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summarize Tissue Variables <a name=\"tissue\"></a>\n",
    "\n",
    "per patient means\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(datadir)\n",
    "#s_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# done\n",
    "# #leiden cell types (fraction)\n",
    "# n_neighbors=30\n",
    "# resolution= 0.6\n",
    "\n",
    "# s_type='NP8-IDC'#'NP9-IDC2'  #\n",
    "# s_subtype = s_type\n",
    "# if s_type == 'NP8-IDC':\n",
    "#     s_sample='20240521_NP-Tissues-ROIs_s'#\n",
    "#     n_markers=19\n",
    "# else:\n",
    "#     n_markers=21# 19#\n",
    "#     s_sample = '20240522_NP9-IDC2-ROIs_s'#\n",
    "# df_lei = pd.read_csv(f'{s_sample}_LeidenClustering_neighbors{n_neighbors}_resolution{resolution}_markers{n_markers}_{s_type}.csv',index_col=0)\n",
    "# df_lei['slide_scene'] = [item.split('_cell')[0] for item in df_lei.index]\n",
    "# df_lei['leidencelltype2'] = df_lei.leidencelltype3.replace({'tumor':'epithelial','endothelial':'stromal','immune':'stromal'})\n",
    "# df_lei['celltype1'] = 'all'\n",
    "# df_lei['countme'] = True\n",
    "\n",
    "# #leiden\n",
    "# for s_celltype in ['leidencelltype2','celltype1']: #'celltype3','celltype',\n",
    "#     for s_cell in df_lei.loc[:,s_celltype].unique():\n",
    "#         df_cell = df_lei.loc[df_lei.loc[:,s_celltype]==s_cell]\n",
    "#         df_prop = viz.prop_positive(df_cell,s_cell='annotated',s_grouper='ROI')\n",
    "#         s_out = f'results_{s_sample}_LeidenClustering_{n_neighbors}_{n_markers}_{resolution}_byPatient_by{s_celltype}_in{s_cell}_{s_type}.csv'\n",
    "#         df_prop.fillna(0).to_csv(s_out)\n",
    "#         print(s_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## compare tumor, border and normal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statannotations.Annotator import Annotator\n",
    "import statsmodels\n",
    "def annotated_stripplot_hue(df,x,y,hue,figsize,hue_order=None,b_correct=True,s=2):\n",
    "    '''\n",
    "    with hue, show pvalues, default FDR corrected \n",
    "    '''\n",
    "    order = df.loc[:,x].unique()\n",
    "    if hue_order ==  None:\n",
    "        hue_order = df.loc[:,hue].unique()\n",
    "    plotting = {\"data\":df,\"x\":x,\"y\":y,\"order\":order,\n",
    "               \"hue\":hue,\"hue_order\":hue_order}\n",
    "    pairs = [((item,hue_order[0]),(item,hue_order[1])) for item in order]\n",
    "\n",
    "    fig,ax = plt.subplots(dpi=200,figsize=figsize)\n",
    "    f,a = plt.subplots()\n",
    "    sns.stripplot(**plotting,dodge=True,ax=ax,s=s,alpha=0.7)\n",
    "    sns.boxplot(**plotting,ax=ax,showmeans=True,medianprops={'visible': False},\n",
    "                                   whiskerprops={'visible': False},meanline=True,showcaps=False,\n",
    "                           meanprops={'color': 'k', 'ls': '-', 'lw': 2},showfliers=False,showbox=False) \n",
    "    annot = Annotator(a,pairs,**plotting,verbose=False)\n",
    "    h, l = ax.get_legend_handles_labels()\n",
    "    ax.legend(h[0:len(hue_order)],l[0:len(hue_order)],bbox_to_anchor=(1,1),title=hue)\n",
    "    annot.configure(test='t-test_ind',text_format=\"star\")\n",
    "    annot.apply_test()\n",
    "    a, test_results = annot.apply_test().annotate()\n",
    "    plt.close(f)\n",
    "    annot = Annotator(ax,pairs,**plotting,verbose=False)\n",
    "    d_pval = {}\n",
    "    for res in test_results:\n",
    "        d_pval.update({res.data.group1[0]:res.data.pvalue})\n",
    "    pvalues = [d_pval[item] for item in order]\n",
    "    reject, corrected, __, __ = statsmodels.stats.multitest.multipletests(pvalues,method='fdr_bh')\n",
    "    formatted_pvalues = [f'p={pvalue:.2}' for pvalue in list(pvalues)]\n",
    "    if b_correct:\n",
    "        formatted_pvalues = [f'FDR={pvalue:.2}' for pvalue in list(corrected)]\n",
    "    annot.set_custom_annotations(formatted_pvalues)\n",
    "    annot.annotate()\n",
    "    return(fig, ax)\n",
    "\n",
    "def name_celltype(str):\n",
    "    str = str.replace('Normal breast luminal','Norm. breast').replace('myoep',\n",
    "        'myoep.').title().replace('Cd','CD').replace('Coli','ColI').replace('Er+','ER+').replace('Fb','Fb.')\n",
    "    return(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "for s_type in ['NP9-IDC2','NP8-IDC',]:\n",
    "    df_file = pd.DataFrame(index=os.listdir())\n",
    "    df_file = df_file[df_file.index.str.contains('results_')]\n",
    "    for hue_order in [('Tumor','Border'),('Tumor','Normal')]:\n",
    "        sns.set_palette('muted')\n",
    "        if hue_order==('Tumor','Normal'):\n",
    "            pal = ('#4878d0','#6acc64')\n",
    "            sns.set_palette(pal)\n",
    "        for s_file in df_file[df_file.index.str.contains(s_type)].index:\n",
    "            df=pd.read_csv(s_file,index_col=0)\n",
    "            s_compartment = s_file.split('_in')[1].split('_')[0]\n",
    "            df_long = df.unstack().reset_index().rename({'level_0':'celltype',0:'Fraction'},axis=1)\n",
    "            df_long['Type'] = [item.split('_')[-1][0:-1] for item in df_long.ROI]\n",
    "            if s_compartment == 'epithelial':\n",
    "                figsize=(5,3)\n",
    "            elif s_compartment == 'stromal':\n",
    "                figsize=(8,3)\n",
    "                if s_type=='NP9-IDC2':\n",
    "                    figsize=(9,3)\n",
    "            else:\n",
    "                 #continue\n",
    "                 figsize=(12,3)\n",
    "            if not s_compartment=='nan':\n",
    "                fig,ax = annotated_stripplot_hue(df_long,'celltype','Fraction','Type',figsize=figsize,\n",
    "                                                 hue_order=hue_order,b_correct=True,s=8)\n",
    "                ax.set_ylabel(f'Fraction in {s_compartment}')\n",
    "                labels = ax.get_xticklabels()\n",
    "                ax.set_xticklabels([name_celltype(item.get_text()) for item in labels])\n",
    "                ax.set_title(s_type)\n",
    "                ax.set_xlabel('')\n",
    "                plt.tight_layout()\n",
    "                fig.savefig(f'figures/stripplot_{s_type}_{s_compartment}.png')\n",
    "        if s_type == 'NP8-IDC':\n",
    "            break\n",
    "    #     break\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # run this on re-registered images to create crops\n",
    "# '/home/groups/graylab_share/Chin_Lab/ChinData/Cyclic_Workflow/cmIF_2019-02-01_NPDCIS2/RegisteredImages/ReRegistered-NP9-IDC2'\n",
    "# from skimage import io, transform\n",
    "# import skimage\n",
    "# import numpy as np\n",
    "# s_file = 'ReRegistered-R1_PCNA.CD8.PD1.CK19_NP9-IDC2-Scene-003_c5_ORG.tif'\n",
    "# img = io.imread(s_file)\n",
    "# a_rescale = skimage.exposure.rescale_intensity(img,in_range=(0,1.5*np.quantile(img,0.9999)))\n",
    "# a_crop = (a_rescale/256).astype(np.uint8)\n",
    "# a_crop = transform.rescale(a_crop, scale=0.5,preserve_range=True,anti_aliasing=True)\n",
    "# io.imsave(f'Cropped/{s_file}',a_crop)"
   ]
  }
 ],
 "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
}