20240528_U54-TMA_cycIF_cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster TMAs\n",
    "\n",
    "**Samples:** \n",
    "- U54 TMA Cyclic\n",
    "\n",
    "**Method**: \n",
    "\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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# # pull data from BCC - done\n",
    "# datadir = f'/home/groups/BCC_Chin_Lab/ChinData/Cyclic_Analysis/cmIF_2021-05-03_PDAC/U54-TMA-9/Jenny/data'\n",
    "# sorted(os.listdir(datadir))"
   ]
  },
  {
   "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/data'\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",
    "1. [load data](#cluster)\n",
    "3. [select normal tissue](#load)\n",
    "5. [single cell Umap](#l7)\n",
    "6. [Leiden celltypes](#l8)\n",
    "7. [annotate leiden celltypes](#lbar)\n",
    "7. [Tissue means](#tissue)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normal Tissues  <a name=\"cluster\"></a>\n",
    "\n",
    "\n",
    "get normal cores from PDAC TMA (done)\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# sorted(os.listdir(f'{datadir}'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# #MIBI int\n",
    "# os.chdir(datadir)\n",
    "# s_sample =  '20230915_U54-TMA_BGSubtractedMeanIntensity.csv'#'20230915_U54-TMA_FilteredMeanIntensity.csv'\n",
    "# df_mi = pd.read_csv(f'{datadir}/{s_sample}',index_col=0)\n",
    "# sorted(df_mi.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_mi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #select\n",
    "\n",
    "# ls_col=['label','fov','centroid-0', 'centroid-1','area', 'eccentricity',\n",
    "#        'Amylase','BCAT', 'CD11c', 'CD14', 'CD163', 'CD20', 'CD206', 'CD209', 'CD3',\n",
    "#        'CD31', 'CD4', 'CD44', 'CD45', 'CD56', 'CD68', 'CD74', 'CD8', 'CD9',\n",
    "#        'CHP', 'COL1A1', 'Calprotectin', 'Chym_Tryp', 'ECAD', 'FAP', 'FN1',\n",
    "#        'FOXP3', 'GATA6', 'HLADR', 'ITGB1', 'KRT17', 'KRT18', 'Ki67', 'MYC',\n",
    "#        'Nuclear', 'PIN1', 'PanKRT', 'Podoplanin', 'SMA', 'Synaptophsyin',\n",
    "#        'TCF7', 'Vimentin',]\n",
    "# df_mi=df_mi.loc[:,ls_col]\n",
    "# df_mi.rename({ 'centroid-1':'DAPIX', 'centroid-0':'DAPIY',},axis=1,inplace=True)"
   ]
  },
  {
   "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({'DAPIX':'DAPI_X','DAPIY':'DAPI_Y'},axis=1,inplace=True)\n",
    "# #df_mi['slide_scene'] = [item.split('_cell')[0] for item in df_mi.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #read in metadata\n",
    "# df_meta = pd.read_csv(f'{codedir}/data/U54-TMA_annotation_update_18Oct2023.csv')\n",
    "# #df_meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# #df_mi['slide_scene'] = df_mi.fov.map(dict(zip(df_meta.MIBI_fov,df_meta.coor_mplexable)))\n",
    "# df_mi['Tissue'] = df_mi.slide_scene.map(dict(zip(df_meta.coor_mplexable,df_meta.Tissue)))\n",
    "# df_mi = df_mi[df_mi.Tissue.str.contains('Normal')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df_mi.Tissue.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %matplotlib inline\n",
    "# fig, ax = plt.subplots()\n",
    "# __ = sns.kdeplot(df_mi.Ecad,ax=ax\n",
    "#            )\n",
    "# ax.axvline(500)\n",
    "# #df_mi = df_mi[df_mi.DAPI11Q>1000]"
   ]
  },
  {
   "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/20240523_Her2B-K154-157_FilteredMeanIntensity.csv'\n",
    "# s_out = f'{codedir}/data/20240528_U54-TMA-Normal-cycIF_BGSubtractedMeanIntensity.csv'\n",
    "# if not os.path.exists(s_out):\n",
    "#     print('saving csv')\n",
    "#     df_mi.to_csv(s_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df_mi['Ecad_Pos'] = df_mi.Ecad>500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ## visualize\n",
    "\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=2,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(f\"{s_scene}\\n{df_pos.loc[:,'Tissue'].iloc[0]}\",y= 0.99)\n",
    "#     return(fig, ax)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_scene = df_mi[df_mi.slide_scene=='U54-TMA-9_sceneE12']\n",
    "fig, ax = plt.subplots()\n",
    "sns.kdeplot(df_scene.CSF1R)\n",
    "ax.axvline(12.5*256)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_int(df_scene,s_color,i_thresh):\n",
    "    #plot\n",
    "    df_scene = df_scene.copy()\n",
    "    s_scene = df_scene.slide_scene.iloc[0]\n",
    "    df_scene.loc[df_scene.loc[:,s_color] < i_thresh,s_color]=0\n",
    "    fig, ax = plt.subplots(figsize=(5,5)) #figsize=(18,12)\n",
    "    #plot cells\n",
    "    sc = ax.scatter(data=df_scene, x='DAPI_X',y='DAPI_Y',c=s_color,s=2,alpha=0.8)\n",
    "    ax.axis('equal')\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    ax.set_title(f'{s_color}')\n",
    "    fig.colorbar(sc)\n",
    "    fig.suptitle(f\"{s_scene}\\n{df_scene.loc[:,'Tissue'].iloc[0]}\",y= 0.99)\n",
    "    return(fig, ax)\n",
    "plot_int(df_scene,'CSF1R',12.5*256)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# #plot ROIS\n",
    "# for s_scene in df_mi.slide_scene.unique():\n",
    "#     df_scene=df_mi[df_mi.slide_scene==s_scene]\n",
    "#     fig, ax = plot_pos(df_scene,['Ecad_Pos'],df_scene) #,'ROI2','ROI3','ROI0','ROI1'\n",
    "#     #break"
   ]
  },
  {
   "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": [
    "s_out = f'{codedir}/data/20240528_U54-TMA-Normal-cycIF_FilteredMeanIntensity.csv'\n",
    "df_mi = pd.read_csv(f'{s_out}',index_col=0) \n",
    "print(len(df_mi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(df_mi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = f'{codedir}/data'\n",
    "os.chdir(datadir)\n",
    "#load normal cores\n",
    "\n",
    "s_out = f'{codedir}/data/20240528_U54-TMA-Normal-cycIF_BGSubtractedMeanIntensity.csv'\n",
    "df_norm  = pd.read_csv(f'{s_out}',index_col=0) \n",
    "# previous BG subtract on CSF1R was too high (it was 14)\n",
    "df_norm.loc[df_mi.index,'CSF1R'] = df_mi.CSF1R\n",
    "print(len(df_norm))\n",
    "d_bg = {'EGFR':4,'CK5':3,'CD3':2,'Vim':2,'CSF1R':12.5}#,\n",
    "#subtract bg - once\n",
    "for s_marker, i_thresh in d_bg.items():\n",
    "    df_norm[s_marker] = np.clip(df_norm.loc[:,s_marker] - i_thresh*256, a_min=0.01, a_max=None) #amin not 0\n",
    "    if s_marker == 'CD3':\n",
    "        df_norm.loc[df_norm.index.str.contains('_sceneG12'),s_marker] = np.clip(df_norm.loc[df_norm.index.str.contains('_sceneG12'),s_marker] - 3*256, a_min=0.01, a_max=None)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sns.kdeplot(df_norm.CSF1R)\n",
    "ax.axvline(1*256)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select features/cells\n",
    "s_sample = s_out.split('data/')[1].split('_Filtered')[0]\n",
    "s_type = 'cycIF'#'Her2B-K154-6S'#'Her2B-K157-6S'#\n",
    "\n",
    "#select markers for cell typing\n",
    "\n",
    "#cyclic\n",
    "ls_col = [# 'FoxP3',#'H3K27', 'H3K4','cPARP', 'gH2AX',\n",
    "       'Ki67', #'LamAC', 'MSH6', 'PCNA', 'RAD51', 'GRNZB', \n",
    "        'CD20','p63', #'pERK', 'pHH3', 'pMYC', 'pRB', 'pRPA', \n",
    "       'CD3', 'CD31', 'CD4', 'CD45', 'CD68', #'CD8', #'ColI', 'ColIV', 'CoxIV',\n",
    "       'PD1', 'CAV1',  'CK19',  'CK7', 'CD90','CK14','CK5',#'CK17','CD44',\n",
    "         'Ecad',   'SYP', 'Vim','CSF1R',# 'PDPN',#'MUC1',# 'TFF1','PDL1','CK8',\n",
    "       'aSMA', 'panCK','EGFR', ]#'pAKT', 'pS6RP',  'area', 'eccentricity',\n",
    "print(len(ls_col))\n",
    "#exclude NAs\n",
    "df_exclude = df_norm[df_norm.loc[:,ls_col].isna().any(axis=1)]\n",
    "adata = sc.AnnData(df_norm.loc[(~df_norm.index.isin(df_exclude.index)),ls_col]) \n",
    "adata.obs['slide_scene'] = adata.obs.index.map(dict(zip(df_norm.index.astype('str'),df_norm.slide_scene)))\n",
    "adata.obs['subtype'] = adata.obs.index.map(dict(zip(df_norm.index.astype('str'),df_norm.Tissue)))\n",
    "\n",
    "if s_type == 'all':\n",
    "    adata = adata#[adata.obs.subtype.isin(['Normal', 'Border', 'Tumor'])]\n",
    "else:\n",
    "    # match to MIBI dataset which lacks normal lung and one placenta core (I12)\n",
    "    adata = adata[(adata.obs.subtype.isin(['Normal_Placenta' ,'Normal_Pancreas', 'Normal_Lymphnode' ,'Normal_Tonsil','Normal_Liver'])) & (~adata.obs.slide_scene.isin(['U54-TMA-9_sceneI12']))]  \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())\n",
    "adata.obs.slide_scene.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'] = adata.obs.index.map(dict(zip(df_norm.index.astype('str'),df_norm.slide_scene)))\n",
    "        adata.obs['subtype'] = adata.obs.index.map(dict(zip(df_norm.index.astype('str'),df_norm.Tissue)))\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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figname = f\"Umap{s_sample}_scene_{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='slide_scene',wspace=.25,save=figname,ax=ax)"
   ]
  },
  {
   "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": {
    "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,title=s_type)\n",
    "\n"
   ]
  },
  {
   "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')\n"
   ]
  },
  {
   "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": [
    "datadir = f'{codedir}/data'\n",
    "os.chdir(f'{datadir}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load ROIS\n",
    "\n",
    "s_out = f'{codedir}/data/20240528_U54-TMA-Normal-cycIF_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",
    "s_type='cycIF'#'Her2B-K154-6S'#\n",
    "\n",
    "if s_type == 'MIBI':\n",
    "    n_markers = 34\n",
    "else:\n",
    "    n_markers = 22 #26 #\n",
    "results_file = f'{s_sample}_{n_neighbors}neighbors_{n_markers}markers_{s_type}_leiden{resolution}.h5ad'\n",
    "results_file='20240528_U54-TMA-Normal-cycIF_BGSubtractedMeanIntensity.csv_s_30neighbors_22markers_cycIF_leiden0.6.h5ad'\n",
    "adata = sc.read_h5ad(results_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_tissue = dict(zip(df_norm.index,df_norm.Tissue))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_norm.Tissue.unique()"
   ]
  },
  {
   "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",
    "df_p['Tissue'] = df_p.index.map(d_tissue)\n",
    "\n",
    "if results_file =='20240528_U54-TMA-Normal-cycIF_BGSubtractedMeanIntensity.csv_s_30neighbors_22markers_cycIF_leiden0.6.h5ad':\n",
    "    print(s_type)\n",
    "    ls_drop = ['24','23','22','14b','13b'] #24 bright CD90, 23 dim, 22 artifact, 13b and 14b atrifact\n",
    "    df_p.loc[(df_p.leiden=='14') & (df_p.Tissue!='Normal_Liver'),'leiden'] = '14b' \n",
    "    df_p.loc[(df_p.leiden=='13') & (df_p.Tissue!='Normal_Liver'),'leiden'] = '13b' \n",
    "    df_p.loc[(df_p.leiden=='8') & (df_p.Tissue=='Normal_Lymphnode'),'leiden'] = '13b' \n",
    "    #df_p.loc[(df_p.leiden=='20') & (df_p.Tissue!='Normal_Tonsil'),'leiden'] = '20b' \n",
    "    df_p.loc[(df_p.leiden=='19') & (df_p.CK5<0.5),'leiden'] = '19t'\n",
    "    df_p.loc[(df_p.leiden=='19t') & (df_p.CD20>1),'leiden'] = '19b'\n",
    "    df_p.loc[(df_p.leiden=='2') & (df_p.CD3<0.5),'leiden'] = '2b'\n",
    "    df_p.loc[(df_p.leiden=='8') & (df_p.CD20>0.5),'leiden'] = '8b'\n",
    "    df_p.loc[(df_p.leiden=='4') & (df_p.CD20>0.5),'leiden'] = '4b'\n",
    "    df_p.loc[(df_p.leiden=='15') & (df_p.CD20>0.5),'leiden'] = '15b'\n",
    "    df_p.loc[(df_p.leiden=='15') & (df_p.index.str.contains('sceneB01')),'leiden'] = '15b'\n",
    "    tum_clust = ['19','18','16','14','9','6','13']\n",
    "    str_clust = ['8']\n",
    "    imm_clust = ['21','20','19b','19t','12','11','10','7','2','2b','8b','5','1','0','4b','15b'] \n",
    "    endo_clust = ['17','3'] #'16'lymphatic\n",
    "    fb_clust = ['4','15']\n",
    "    d_named = {'21':'TFH cell','20':'B cell','19b':'B cell','19':'Tonsil Crypt','19t':'T cell',\n",
    "               '18':'Ductal','17':'aSMA+ blood vessel','16':'acinar','15':'p63+ str.',\n",
    "               '14':'Liver lobule','13':'Liver lobule','12':'CD4 T cell','11':'CD68+ mac.',\n",
    "               '10':'CD4 T cell','7':'CD4 T cell','2':'CD4 T cell','2b':'B cell','15b':'B cell',\n",
    "               '9':'acinar','8b':'B cell','8':'Quies. str.','6':'acinar','5':'Ki67+ B cell',\n",
    "               '4':'Vim+ FB','4b':'B cell','3':'endothelial','1':'B cell','0':'CD4 T cell'\n",
    "              }\n",
    "\n",
    "else:\n",
    "    print('new one')\n",
    "    ls_drop=[]\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 +endo_clust+str_clust +ls_drop+fb_clust)\n",
    "#set(tum_clust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, ax =plt.subplots()\n",
    "# sns.kdeplot(df_p.CD3,ax=ax)\n",
    "# ax.axvline(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clustermap\n",
    "if s_type=='Her2B-K154-6S':\n",
    "    ls_drop_col = []\n",
    "else:\n",
    "    ls_drop_col = ['Tissue','SYP']\n",
    "b_annot = True#True#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,8),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.fig.suptitle(s_type)\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,4),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": {},
   "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": [
    "\n",
    "s_clust = '8'\n",
    "df_p['slide_scene'] = df_p.index.map(dict(zip(df_norm.index.astype('str'),df_norm.slide_scene)))#[item.split('_cell')[0] for item in df_p.index]\n",
    "\n",
    "ls_sene = df_p.groupby(['leiden','slide_scene']).count().loc[s_clust].iloc[:,1].sort_values(ascending=False).index[0:5]\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": [
    "#plot all groups spatially - leiden \n",
    "ls_clust = [s_clust,'15b']#df_pos.columns.tolist()#['8','11','12','7','5','10','1','9']#'0','2',\n",
    "for s_slide in ls_sene:#sorted(set(df_p.slide_scene)):\n",
    "    fig,ax = plt.subplots(figsize=(5,4.5),dpi=200)\n",
    "    #plot negative cells\n",
    "    df_scene = df_data[df_data.slide_scene.str.contains(s_slide)].copy()\n",
    "    df_scene.index = df_scene.index.astype('str')\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",
    "            print(len(ls_index))\n",
    "            ax.scatter(data=df_scene.loc[ls_index],x='DAPI_X',y='DAPI_Y',label=f'{s_color}',s=0.1,color=colors[idxs])\n",
    "        #break\n",
    "    \n",
    "    ax.set_title(f\"{s_slide} {df_norm.loc[df_norm.slide_scene==s_slide,'Tissue'].iloc[0]}\", 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','aSMA','Ecad','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": [
    "#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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_p['annotated'] = df_p.leiden.map(d_named).fillna('_Artifact')\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 = 'U54-TMA-9_sceneG06'\n",
    "    fig,ax = plt.subplots(figsize=(8,8),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.5,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.5,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=8,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  os.path.exists(s_out):\n",
    "    print(f'saving csv {s_out}')\n",
    "    df_out.to_csv(s_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#s_out = '20240523_Her2B-K154-157_s_LeidenClustering_neighbors30_resolution0.6_markers25_Her2B-K157-6S.csv'\n",
    "#s_out = '20240523_Her2B-K154-157_s_LeidenClustering_neighbors30_resolution0.6_markers24_Her2B-K154-6S.csv'\n",
    "df_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_out = '20240523_Her2B-K154-157_s_LeidenClustering_neighbors30_resolution0.6_markers25_Her2B-K157-6S.csv'\n",
    "#s_out = '20240523_Her2B-K154-157_s_LeidenClustering_neighbors30_resolution0.6_markers24_Her2B-K154-6S.csv'\n",
    "df_lei['matched']=df_lei.leidencelltype4.copy()\n",
    "df_lei.loc[df_lei.annotated.str.contains('B cell'),'matched'] = 'B cell'\n",
    "df_lei.loc[df_lei.annotated.str.contains('T cell'),'matched'] = 'T cell'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample ='20240528_U54-TMA-Normal-cycIF_s'\n",
    "#leiden cell types (fraction)\n",
    "n_neighbors=30\n",
    "resolution= 0.6\n",
    "\n",
    "s_type= 'cycIF'#'Her2B-K157-6S'#\n",
    "s_subtype = s_type\n",
    "if s_type == 'Her2B-K157-6S':\n",
    "    n_markers=25\n",
    "else:\n",
    "    n_markers=22\n",
    "\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",
    "df_lei['matched']=df_lei.leidencelltype4.copy()\n",
    "df_lei.loc[df_lei.annotated.str.contains('B cell'),'matched'] = 'B cell'\n",
    "df_lei.loc[df_lei.annotated.str.contains('T cell'),'matched'] = 'T cell'\n",
    "df_lei.loc[df_lei.annotated == 'Ductal','matched'] = 'Ductal'\n",
    "#leiden\n",
    "#leiden\n",
    "for s_celltype in ['celltype1']: #'celltype3','celltype','leidencelltype2',\n",
    "    for s_cell in df_lei.loc[:,s_celltype].unique():\n",
    "        for s_grouper in ['slide_scene','celltype1']:\n",
    "            df_cell = df_lei.loc[df_lei.loc[:,s_celltype]==s_cell]\n",
    "            df_prop = viz.prop_positive(df_cell,s_cell='matched',s_grouper=s_grouper)\n",
    "            s_out = f'results_{s_sample}_LeidenClustering_{n_neighbors}_{n_markers}_{resolution}_by{s_grouper}_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": {},
   "outputs": [],
   "source": [
    "#d_tissue = dict(zip(df_lei.slide_scene, df_lei.Tissue))\n",
    "d_tissue= {'U54-TMA-9_sceneA01': 'Normal_Placenta',\n",
    " 'U54-TMA-9_sceneC07': 'Normal_Liver',\n",
    " 'U54-TMA-9_sceneE12': 'Normal_Liver',\n",
    " 'U54-TMA-9_sceneG12': 'Normal_Pancreas',\n",
    " 'U54-TMA-9_sceneH02': 'Normal_Pancreas',\n",
    " 'U54-TMA-9_sceneB01': 'Normal_Lymphnode',\n",
    " 'U54-TMA-9_sceneG06': 'Normal_Lymphnode',\n",
    " 'U54-TMA-9_sceneD03': 'Normal_Tonsil',\n",
    " 'U54-TMA-9_sceneF08': 'Normal_Tonsil'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = mpl.cm.tab10.colors + mpl.cm.Accent.colors \n",
    "for s_compartment in ['inall']:\n",
    "    df_file = pd.DataFrame(index=os.listdir())\n",
    "    df_file = df_file[(df_file.index.str.contains('results_20240528_U54-TMA-Normal')) & (df_file.index.str.contains('byslide_scene'))]\n",
    "    df_all = pd.DataFrame()\n",
    "    ls_index = sorted(df_file[df_file.index.str.contains(s_compartment)].index)\n",
    "    df1=pd.read_csv(ls_index[0],index_col=0)\n",
    "    df1.index = [item.split('_')[-1] for item in df1.index]\n",
    "    df2=pd.read_csv(ls_index[1],index_col=0)\n",
    "    df2.index = [item.split('_')[-1] for item in df2.index]\n",
    "    s_1 = ls_index[0].split('_')[-1].split('.')[0]\n",
    "    print(s_1)\n",
    "    s_2 = ls_index[1].split('_')[-1].split('.')[0]\n",
    "    fig, ax = plt.subplots(3,3,figsize=(8,10),sharex=True,sharey=True,dpi=300)\n",
    "    ax=ax.ravel()\n",
    "    ls_r = []\n",
    "    for idx, s_index in enumerate(df1.index):\n",
    "        x = df1.loc[s_index]\n",
    "        y = df2.loc[s_index]\n",
    "        r, p = scipy.stats.pearsonr(x.values, y.values)\n",
    "        ls_r.append(r)\n",
    "        df = pd.concat([pd.DataFrame(x),pd.DataFrame(y)],axis=1)\n",
    "        s_tissue = d_tissue['U54-TMA-9_'+df.columns[0]].replace(\"_\",\" \")\n",
    "        s_scene = df.columns[0].replace(\"scene\",\"\")#.replace(\"00\",\" \")\n",
    "        df.columns = [s_1,s_2]\n",
    "        df.rename({'AF':'_AF'},inplace=True)\n",
    "        #plot\n",
    "        title = f'{s_scene.title()}: {s_tissue.title()}\\nr={r:.3} p={p:.2}'\n",
    "        df.T.plot(kind='bar',stacked=True,ax=ax[idx],legend=None,title=title,color=colors)\n",
    "        if idx < 6:\n",
    "            ax[idx].tick_params(\n",
    "                axis='x',          # changes apply to the x-axis\n",
    "                which='both',      # both major and minor ticks are affected\n",
    "                bottom=False,      # ticks along the bottom edge are off\n",
    "                top=False,         # ticks along the top edge are off\n",
    "                labelbottom=False) # labels along the bottom edge are off\n",
    "    ax[idx].legend(bbox_to_anchor=(1,1),ncol=2,title='Cell Type')\n",
    "\n",
    "    ax[0].set_ylabel(f'Fraction in {s_compartment.split(\"in\")[1]}')\n",
    "    ax[3].set_ylabel(f'Fraction in {s_compartment.split(\"in\")[1]}')\n",
    "    ax[6].set_ylabel(f'Fraction in {s_compartment.split(\"in\")[1]}')\n",
    "    #ax[3].set_xticklabels(['154','157'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.array(ls_r).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.array(ls_r).std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#all\n",
    "colors = mpl.cm.tab10.colors + mpl.cm.Accent.colors \n",
    "for s_compartment in ['inall']:\n",
    "    df_file = pd.DataFrame(index=os.listdir())\n",
    "    df_file = df_file[(df_file.index.str.contains('results_20240528_U54-TMA-Normal')) & (df_file.index.str.contains('_bycelltype1_bycelltype1'))]\n",
    "    df_all = pd.DataFrame()\n",
    "    ls_index = sorted(df_file[df_file.index.str.contains(s_compartment)].index)\n",
    "    df1=pd.read_csv(ls_index[0],index_col=0)\n",
    "    df1.index = [item.split('_')[-1] for item in df1.index]\n",
    "    df2=pd.read_csv(ls_index[1],index_col=0)\n",
    "    df2.index = [item.split('_')[-1] for item in df2.index]\n",
    "    s_1 = ls_index[0].split('_')[-1].split('.')[0]\n",
    "    print(s_1)\n",
    "    s_2 = ls_index[1].split('_')[-1].split('.')[0]\n",
    "    fig, ax = plt.subplots(figsize=(3,4),sharex=True,sharey=True,dpi=300)\n",
    "    ls_r = []\n",
    "    for idx, s_index in enumerate(df1.index):\n",
    "        x = df1.loc[s_index]\n",
    "        y = df2.loc[s_index]\n",
    "        r, p = scipy.stats.pearsonr(x.values, y.values)\n",
    "        ls_r.append(r)\n",
    "        df = pd.concat([pd.DataFrame(x),pd.DataFrame(y)],axis=1)\n",
    "        s_tissue = '' #d_tissue['U54-TMA-9_'+df.columns[0]].replace(\"_\",\" \")\n",
    "        s_scene = df.columns[0].replace(\"scene\",\"\")#.replace(\"00\",\" \")\n",
    "        df.columns = [s_1,s_2]\n",
    "        df.rename({'AF':'_AF'},inplace=True)\n",
    "        #plot\n",
    "        title = f'Adjacent Section Correlation\\nr={r:.2} p={p:.2}'\n",
    "        df.T.plot(kind='bar',stacked=True,ax=ax,legend=None,title=title,color=colors)\n",
    "        ax.set_title(title,fontsize=16)\n",
    "    ax.legend(bbox_to_anchor=(1,1),ncol=2,title='Cell Type')"
   ]
  },
  {
   "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
}