20220410_MIBI_cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster MIBI Samples\n",
    "\n",
    "**Samples:** \n",
    "- MIBI FOV Size: 800 μm Frame Size: 2048 (0.39 μm per pixel)\n",
    "\n",
    "**Method**: Cluster Umap with Leiden algorithm \n",
    "\n",
    "Compare to images\n",
    "\n",
    "compare to hitological subtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load libraries\n",
    "import pandas as pd\n",
    "import re\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 matplotlib.pyplot as plt\n",
    "mpl.rc('figure', max_open_warning = 0)\n",
    "codedir = os.getcwd()\n",
    "import phenograph\n",
    "os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data')\n",
    "from mplex_image import visualize as viz, process, preprocess, normalize, mics, mpimage\n",
    "np.random.seed(1027)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "codedir = '/home/groups/graylab_share/OMERO.rdsStore/engje/Data/Keren et al_'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change to correct directory\n",
    "datadir =  f'{codedir}'\n",
    "s_date = '20220410'\n",
    "os.chdir(codedir)\n",
    "if not os.path.exists(s_date):\n",
    "    os.mkdir(s_date)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.listdir()\n",
    "#os.listdir('results')\n",
    "#os.listdir('intermediate_data/protein_positivity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of contents <a name=\"contents\"></a>\n",
    "2. [Load annotation](#loadold)\n",
    "2. [Load intensity](#load)\n",
    "4. [Normalize](#leiden)\n",
    "5. [Umap](#l7)\n",
    "6. [Leiden](#l8)\n",
    "[annotate](#lbar)\n",
    "5. [Sub-cluster](#kmeans)\n",
    "8. [Visualize](#viz)\n",
    "9. [subtype](#subtype)\n",
    "9. [survival](#surv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Annotation <a name=\"loadold\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#positve\n",
    "df_name = pd.read_csv(f'{datadir}/proteins_by_frame.csv',index_col=0)\n",
    "df_a = pd.read_csv(f'{datadir}/clinical_data.csv',index_col=0)\n",
    "#load annotation\n",
    "df_a['Group'] = pd.qcut(x=df_a.loc[:,'Survival_time'],q=2,labels=['short','long'])\n",
    "d_a = dict(zip(df_a.index.astype('str'),df_a.Group))\n",
    "d_name = dict(zip(df_name.index.astype('str'),df_name.Biomarker))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Mean Intensity <a name=\"load\"></a>\n",
    "\n",
    "load data\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi_full = pd.read_csv(f'{codedir}/Segmentation/features_MIBI-TMA_MeanIntensity_Centroid_Shape.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi_full.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set([item.split('_')[0] + '_cytoplasm' for item in sorted(df_mi_full.columns)])\n",
    "ls_marker = [#'Au_cytoplasm', 'Background_cytoplasm',\n",
    " 'Beta catenin_cytoplasm','Beta catenin_nuclei',\n",
    " 'CD11b_cytoplasm', 'CD11c_cytoplasm', 'CD138_cytoplasm', 'CD16_cytoplasm', 'CD209_cytoplasm', 'CD20_cytoplasm',\n",
    " 'CD31_cytoplasm', 'CD3_cytoplasm', 'CD45RO_cytoplasm', 'CD45_cytoplasm', 'CD4_cytoplasm', 'CD56_cytoplasm',\n",
    " 'CD63_cytoplasm', 'CD68_cytoplasm', 'CD8_cytoplasm', #'Ca_cytoplasm','Fe_cytoplasm',\n",
    " 'EGFR_cytoplasm', 'FoxP3_nuclei', 'H3K27me3_nuclei', 'H3K9ac_nuclei',\n",
    " 'HLA-Class-1_cytoplasm', 'HLA-DR_cytoplasm', 'IDO_cytoplasm', 'Keratin17_cytoplasm', 'Keratin6_cytoplasm',\n",
    " 'Ki67_nuclei', 'Lag3_cytoplasm', 'MPO_cytoplasm',# 'Na_cytoplasm',\n",
    " 'PD-L1_cytoplasm', 'PD1_cytoplasm',# 'P_cytoplasm',\n",
    " 'Pan-Keratin_cytoplasm', 'SMA_cytoplasm',# 'Si_cytoplasm', 'Ta_cytoplasm',\n",
    " 'Vimentin_cytoplasm', 'area_nuclei', 'centroid-0_nuclei', 'centroid-1_nuclei', 'eccentricity_nuclei',# 'label_cytoplasm',\n",
    " 'p53_nuclei', 'phospho-S6_cytoplasm']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_replace = {'centroid-0_nuclei':'DAPI_Y',\n",
    "       'centroid-1_nuclei':'DAPI_X',\n",
    "            'Beta catenin_cytoplasm':'Beta-catenin-c',\n",
    "           'Beta catenin_nuclei':'Beta-catenin-n',\n",
    "            'CD11b_cytoplasm': 'CD11b',\n",
    " 'CD11c_cytoplasm': 'CD11c',\n",
    " 'CD138_cytoplasm': 'CD138',\n",
    " 'CD16_cytoplasm': 'CD16',\n",
    " 'CD209_cytoplasm': 'CD209',\n",
    " 'CD20_cytoplasm': 'CD20',\n",
    " 'CD31_cytoplasm': 'CD31',\n",
    " 'CD3_cytoplasm': 'CD3',\n",
    " 'CD45RO_cytoplasm': 'CD45RO',\n",
    " 'CD45_cytoplasm': 'CD45',\n",
    " 'CD4_cytoplasm': 'CD4',\n",
    " 'CD56_cytoplasm': 'CD56',\n",
    " 'CD63_cytoplasm': 'CD63',\n",
    " 'CD68_cytoplasm': 'CD68',\n",
    " 'CD8_cytoplasm': 'CD8',\n",
    " 'EGFR_cytoplasm': 'EGFR',\n",
    " 'FoxP3_nuclei': 'FoxP3',\n",
    " 'H3K27me3_nuclei': 'H3K27me3',\n",
    " 'H3K9ac_nuclei': 'H3K9ac',\n",
    " 'HLA-Class-1_cytoplasm': 'HLA-Class-1',\n",
    " 'HLA-DR_cytoplasm': 'HLA-DR',\n",
    " 'IDO_cytoplasm': 'IDO',\n",
    " 'Keratin17_cytoplasm': 'Keratin17',\n",
    " 'Keratin6_cytoplasm': 'Keratin6',\n",
    " 'Ki67_nuclei': 'Ki67',\n",
    " 'Lag3_cytoplasm': 'Lag3',\n",
    " 'MPO_cytoplasm': 'MPO',\n",
    " 'PD-L1_cytoplasm': 'PD-L1',\n",
    " 'PD1_cytoplasm': 'PD1',\n",
    " 'Pan-Keratin_cytoplasm': 'Pan-Keratin',\n",
    " 'SMA_cytoplasm': 'SMA',\n",
    " 'Vimentin_cytoplasm': 'Vimentin',\n",
    " 'area_nuclei': 'area',\n",
    " 'eccentricity_nuclei': 'eccentricity',\n",
    " 'p53_nuclei': 'p53',\n",
    " 'phospho-S6_cytoplasm': 'phospho-S6'}\n",
    "df_mi = df_mi_full.loc[:,ls_marker].rename(d_replace,axis=1)\n",
    "             "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''['Au',\n",
    " 'Background',\n",
    " 'Beta catenin',#good\n",
    " 'CD11b', #good\n",
    " 'CD11c',#good\n",
    " 'CD138', #? CD45 negative\n",
    " 'CD16', #good\n",
    " 'CD20', #good\n",
    " 'CD209',#? rare DC-Sign\n",
    " 'CD3',#good\n",
    " 'CD31', #okay\n",
    " 'CD4', #good\n",
    " 'CD45',#good\n",
    " 'CD45RO',#same as CD45\n",
    " 'CD56',#rare\n",
    " 'CD63',#good\n",
    " 'CD68',#good\n",
    " 'CD8', #good\n",
    " 'Ca',\n",
    " 'Celltype',\n",
    " 'Column',\n",
    " 'EGFR', #dim, but tumor specific\n",
    " 'Fe',\n",
    " 'FoxP3', #rare\n",
    " 'H3K27me3',#everywhere\n",
    " 'H3K9ac',\n",
    " 'HLA-DR',#good\n",
    " 'HLA_Class_1',#good\n",
    " 'IDO',#probably okay.\n",
    " 'Keratin17',#good\n",
    " 'Keratin6',#good\n",
    " 'Ki67',#good\n",
    " 'Lag3',#rare\n",
    " 'MPO',#rare\n",
    " 'Na',\n",
    " 'P',\n",
    " 'PD-L1',#okay\n",
    " 'PD1',#good rare\n",
    " 'Pan-Keratin',#good\n",
    " 'Row',\n",
    " 'SMA',#brigt\n",
    " 'Sample',\n",
    " 'Si',\n",
    " 'Ta',\n",
    " 'Vimentin',#good\n",
    " 'dsDNA',\n",
    " 'p53',#good\n",
    " 'phospho-S6'#okay\n",
    "]'''\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_drop = ['Au', 'Background', 'Ca',  'Fe', 'Na', 'P','Si', 'Ta'] #'Celltype', 'Column','Row','Sample' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_mi.columns = [item.split('_')[0] for item in df_mi.columns]\n",
    "\n",
    "# Failed':{'CC3-647','CD4','CD3','F4-80','MHCII','ELN','NFkB','CC3-488', 'panCK','Arg1',}\n",
    "df_filter_mi = df_mi.loc[:,~df_mi.columns.isin(ls_drop)]\n",
    "#df_mi.rename({'nuclei_area':'area','nuclei_eccentricity':'eccentricity','slide_scene':'slidescene'},axis=1,inplace=True)\n",
    "#df_xy.rename({'DAPI_X':'DAPI_X_old','DAPI_Y':'DAPI_Y_old','slide_scene':'slide_scene_old'},axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_mi.columns = [item.split('_')[0] for item in df_mi.columns]\n",
    "df_filter_mi.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#s_date\n",
    "s_sample='MIBI'\n",
    "s_date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_out = f'{codedir}/{s_date}_{s_sample}_FilteredMeanIntensity.csv' #1027 old, 1029 updated\n",
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_filter_mi.to_csv(s_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalize <a name=\"leiden\"></a>\n",
    "\n",
    "background subtraction\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "os.chdir(codedir)\n",
    "s_date = '20211027' #326562 rows\n",
    "s_date = '20211029' #748728  rows\n",
    "s_date = '20211109' #722760 rows\n",
    "s_sample = f'{s_date}_RS-mTMA-5'\n",
    "df_mi = pd.read_csv(f'{s_sample}_FilteredMeanIntensity.csv',index_col=0)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subtract background\n",
    "'''\n",
    "d_bg = {'CD31':8,#11,\n",
    "        'CD45':3, 'CD8':10, 'ColIV':12, 'FoxP3':18,#12,\n",
    " 'CD103':10, 'CD11b':1.5, 'CD11c':6\n",
    "}\n",
    "df_bg = df_mi.copy()\n",
    "for s_marker, i_thresh in d_bg.items():\n",
    "    df_bg[s_marker] = np.clip(df_mi.loc[:,s_marker] - i_thresh*256, a_min=0, a_max=None)\n",
    "\n",
    "s_out = f'{codedir}/{s_sample}_BGSubtractedMeanIntensity.csv'\n",
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_bg.to_csv(s_out)\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Umap <a name=\"l7\"></a>\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load annotation\n",
    "df_a['Group'] = pd.qcut(x=df_a.loc[:,'Survival_time'],q=2,labels=['short','long'])\n",
    "d_a = dict(zip(df_a.index.astype('str'),df_a.Group))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "s_sample = '20220316_MIBI'\n",
    "s_date = s_sample.split('_')[0]\n",
    "#22 features cd11b monocytes, granulocytes, macrophages, and natural killer cells\n",
    "#cd11c high levels on most human dendritic cells, but also on monocytes, macrophages, neutrophils, and some B cells \n",
    "lls_col = [     ['Beta-catenin-c', 'Beta-catenin-n', 'CD11b', 'CD11c', 'CD138', 'CD16', 'CD20', 'CD209',\n",
    "       'CD3', 'CD31', 'CD4', 'CD45', 'CD45RO', 'CD56', 'CD63', 'CD68', 'CD8','HLA-Class-1',\n",
    "        'EGFR', 'FoxP3', 'HLA-DR',  'IDO', # 'H3K27me3', 'H3K9ac', 'phospho-S6','HLA',\n",
    "       'Keratin17', 'Keratin6', 'Ki67', 'Lag3', 'MPO', 'p53', 'Pan-Keratin','area','eccentricity',\n",
    "       'PD-L1', 'PD1', 'SMA', 'Vimentin'], #'dsDNA', oops not in new features\n",
    "]\n",
    "\n",
    "for s_type in ['Filtered']:\n",
    "    df_mi = pd.read_csv(f'{s_sample}_{s_type}MeanIntensity.csv',index_col=0)\n",
    "    break\n",
    "    for ls_col in lls_col:\n",
    "        print(len(ls_col))\n",
    "        adata = sc.AnnData(df_mi.loc[:,ls_col])\n",
    "        #raw\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)}_{s_type}.png')\n",
    "        plt.tight_layout()\n",
    "        sc.pl.pca_variance_ratio(adata,n_pcs=32, log=True)\n",
    "        adata.obs['slide'] = [item.split('_')[0] for item in adata.obs.index]\n",
    "        adata.obs['slide'] = adata.obs.slide.replace({'t1':'1'})\n",
    "        sc.pp.scale(adata, zero_center=False, max_value=20)\n",
    "        s_sample = s_sample + '_s'\n",
    "\n",
    "        #load annotation\n",
    "        # mapping\n",
    "        d_a = dict(zip(df_a.index.astype('str'),df_a.Group))\n",
    "        adata.obs['Group'] = adata.obs.slide.map(d_a).fillna('unknown')\n",
    "\n",
    "        #umap\n",
    "        for n_neighbors in [30]: #,30\n",
    "            results_file = f'{s_date}/{s_sample}_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.h5ad'\n",
    "            if os.path.exists(results_file):\n",
    "                print('loading umap')\n",
    "                adata = sc.read_h5ad(results_file)\n",
    "            else:\n",
    "                print('calc umap')\n",
    "                # calculate neighbors \n",
    "                sc.pp.neighbors(adata, n_neighbors=n_neighbors) #, method='rapids'\n",
    "                sc.tl.umap(adata)\n",
    "            #save results\n",
    "            if not os.path.exists(results_file):\n",
    "                adata.write(results_file)\n",
    "            #color by markers\n",
    "            \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.5',ncols=6,save=figname)\n",
    "            #color by slide\n",
    "            figname = f\"Umap_slide_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.png\"\n",
    "            title=figname.split('.png')[0].replace('_',' ')\n",
    "            fig,ax = plt.subplots(figsize=(3,3), dpi=200)\n",
    "            sc.pl.umap(adata, color='slide',wspace=.25,save=figname,ax=ax,title=title)\n",
    "            \n",
    "            #color by group\n",
    "            \n",
    "            figname = f\"Umap_Group_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.png\"\n",
    "            title=figname.split('.png')[0].replace('_',' ')\n",
    "            fig,ax = plt.subplots(figsize=(3,3), dpi=200)\n",
    "            sc.pl.umap(adata, color='Group',wspace=.25,save=figname,ax=ax,title=title)\n",
    "            \n",
    "            break\n",
    "        break\n",
    "    break\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs.slide.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#screen all\n",
    "for s_slide in sorted(adata.obs.slide.unique()):\n",
    "    adata.obs['test'] = (adata.obs.slide == s_slide).replace({True:1,False:0})\n",
    "    sc.pl.umap(adata, color='test',title=s_slide)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# weird ones (fixed)\n",
    "'''\n",
    "s_slide = 'RS-mTMA-5_sceneE03' #corner extra\n",
    "s_slide = 'RS-mTMA-5_sceneH11' # circel imaging artifact (pmyc)\n",
    "s_slide = 'RS-mTMA-5_sceneI02' #cd31 bg\n",
    "#Still unique\n",
    "s_slide = 'RS-mTMA-5_sceneA05' #high foxp3 bg, some float\n",
    "s_slide = 'RS-mTMA-5_sceneE07' #folded tissue\n",
    "s_slide = 'RS-mTMA-5_sceneH06' #bright floating tissue\n",
    "s_slide = 'RS-mTMA-5_sceneI01' #bright floating tissue\n",
    "s_slide = 'RS-mTMA-5_sceneI11' #epcam, ecad, CK5 neg (no tumor)\n",
    "s_slide = 'RS-mTMA-5_sceneI07' #dried out bottom? pmyc weird\n",
    "#new ouliers\n",
    "#A05, E06, E07, F03, G09, H06, I01, I07, I11\n",
    "#other tissues\n",
    "# A11 spleen B10 spleen?\n",
    "#B6 lung #I10 lung\n",
    "#d4 normal pad, d9 adj #d12 fat pad\n",
    "# E6, E12 lymph node I12 lymph node\n",
    "#for s_slide in adata.obs.slide.unique():\n",
    "adata.obs['test'] = (adata.obs.slide == s_slide).replace({True:1,False:0})\n",
    "sc.pl.umap(adata, color='test',title=s_slide)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster Leiden <a name=\"l8\"></a>\n",
    "\n",
    "main cell types\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '20220316_MIBI_s'\n",
    "s_type = 'Filtered'\n",
    "s_date = s_sample.split('_')[0]\n",
    "n_neighbors = 30\n",
    "n_markers = 35\n",
    "\n",
    "results_file = f'{s_date}/{s_sample}_{n_neighbors}neighbors_{n_markers}markers_{s_type}.h5ad'\n",
    "adata = sc.read_h5ad(results_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#leiden\n",
    "for resolution in [0.6,0.4, 0.5, 0.6]: #0.4,\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",
    "        sc.tl.leiden(adata,resolution=resolution)\n",
    "        adata.write(results_file)\n",
    "    else:\n",
    "        print('loading')\n",
    "        adata = sc.read_h5ad(results_file)\n",
    "    fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "    figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}.png'\n",
    "    sc.pl.umap(adata, color='leiden',ax=ax,title=figname.split('.png')[0].replace('_',' '),wspace=.25,save=figname,)\n",
    "    #fig.savefig(f'./figures/{figname}')\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annotate Leiden <a name=\"lbar\"></a>\n",
    "\n",
    "main cell types\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '20220316_MIBI_s'\n",
    "s_type = 'Filtered'\n",
    "s_date = s_sample.split('_')[0]\n",
    "n_neighbors = 30\n",
    "n_markers = 35\n",
    "resolution = 0.6 #0.7\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": [
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}.png'\n",
    "sc.pl.umap(adata, color='leiden',ax=ax,title=figname.split('.png')[0].replace('_',' '),wspace=.25,save=figname,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_col = ['CD31','Beta-catenin-c','Pan-Keratin','Vimentin','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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_p = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[df_p.CD31>.2].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='CD31')\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[df_p.CD45>1].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='CD45')\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[(df_p.loc[:,'Beta-catenin-c']>.4) | (df_p.loc[:,'Pan-Keratin']>1)].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='Beta-catenin-c and Pan-Keratin')\n",
    "\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[(df_p.Vimentin>6)].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='Vimentin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "df_p.columns = [item.replace(' ','_').replace('-','_') for item in df_p.columns]\n",
    "df_p['leiden'] = adata.obs['leiden'].astype('str')\n",
    "#CD45 above 300\n",
    "#CK5 above 2048 (768)\n",
    "#ecad/epcam above 2000 (768)\n",
    "# CD68 above 1536\n",
    "#area above 318\n",
    "#all 30 neigh, 7 pcs\n",
    "if results_file == '20220316_MIBI_s_30neighbors_35markers_Filtered_leiden0.4.h5ad':\n",
    "    tum_clust = ['5', '13','12', '9', '7', '3', '1', '6','10']\n",
    "    str_clust = ['17']\n",
    "    imm_clust = ['4', '2', '14','8','15','16'] \n",
    "    endo_clust = ['11']\n",
    "    fb_clust= ['0']\n",
    "if results_file == '20220316_MIBI_s_30neighbors_35markers_Filtered_leiden0.7.h5ad':\n",
    "    tum_clust = ['8', '18','15', '2', '4', '9', '17','21','12', '6', '3', '22','13']\n",
    "    str_clust = ['16','11']\n",
    "    imm_clust = ['5', '1','23', '10', '19' ,'7', '20'] \n",
    "    endo_clust = ['14']\n",
    "    fb_clust= ['0']\n",
    "    ls_drop= []\n",
    "if results_file == '20220316_MIBI_s_30neighbors_35markers_Filtered_leiden0.6.h5ad':\n",
    "    print('goodone')\n",
    "    #split cd4 cd8\n",
    "    df_p.loc[((df_p.leiden=='2') & (df_p.CD8 > 1.5*df_p.CD4)),'leiden'] = '2b'\n",
    "    #split cd4 cd20 & cd11c\n",
    "    df_p.loc[((df_p.leiden=='3') & (df_p.CD20 > 2.5*df_p.CD4)),'leiden'] = '3b'\n",
    "    df_p.loc[((df_p.leiden=='3') & (df_p.CD11c > df_p.CD4)),'leiden'] = '3c'\n",
    "    #split 14\n",
    "    df_p.loc[((df_p.leiden=='14') & (df_p.Pan_Keratin<.5)),'leiden'] = '14b'\n",
    "    tum_clust = ['14', '7', '15','9', '6', '5', '19', '11','0','8','13']\n",
    "    str_clust = ['21','4']\n",
    "    imm_clust = ['20','3', '2', '17','18','10','16','2b','14b','3b','3c'] \n",
    "    endo_clust = ['12']\n",
    "    fb_clust= ['1']\n",
    "    ls_drop= ['21']\n",
    "\n",
    "\n",
    "#df_p.groupby('leiden').mean().CD45.sort_values(ascending=False).index\n",
    "#set(df_p.groupby('leiden').mean().CD31.sort_values(ascending=False).index) - set(tum_clust + imm_clust + endo_clust)\n",
    "#set(tum_clust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_p[df_p.leiden=='3'].CD20.median()\n",
    "#df_p[df_p.leiden=='3'].CD4.median()\n",
    "#df_p['slide'] = [item.split('_')[0] for item in df_p.index]\n",
    "#s_clust = '3' #11\n",
    "#df_p[df_p.leiden==s_clust].groupby('slide').count().sort_values(by='CD31',ascending=False)[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "d_named = {'3':'CD4 T cell',\n",
    "'3b':'CD20 B cell',\n",
    " '3c':'Dendritic cell',\n",
    " '2b':'CD8 T cell',\n",
    " '2':'CD4 T cell',\n",
    " '17':'FoxP3 Treg',\n",
    " '7':'Basal t.',\n",
    " '15':'PD-L1+ Basal t.',\n",
    " '14':'IDO+ Basal t.',\n",
    " '14b':'Dendritic cell',\n",
    " '18':'NK cell',\n",
    " '12':'endothelial',\n",
    " '20':'CD209+ imm.', #DC-SIGN\n",
    " # artifact'21':'quiesc str',\n",
    " '9':'Prolif. t.',\n",
    " '6':'CK lo. t.',\n",
    " '5':'Myoepithelial',\n",
    " '19':'EGFR+ t.',\n",
    " '11':'CD63+ t.',\n",
    " '13':'CK lo. t.',\n",
    " '1':'FB',\n",
    " '4':'Quies. str.',\n",
    " '0':'Luminal t.',\n",
    " '8':'Basal t.',\n",
    " '10':'Macrophage',\n",
    " '16':'neutrophil'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#clustermap\n",
    "#df_p = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "b_annot=False\n",
    "b_annot = True\n",
    "\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[~df_p.leiden.isin(ls_drop)].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=(11,7),cmap='viridis',row_colors=row_colors,#method='single',\n",
    "                   vmin=-2,vmax=2)\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)\n",
    "else:\n",
    "    g = sns.clustermap(df_plot,z_score=1,figsize=(9,7),cmap='viridis',#row_colors=row_colors,#method='single',\n",
    "                   vmin=-2,vmax=2)\n",
    "g.savefig(f'{codedir}/{s_date}/{s_sample}_clustermap_leiden_{resolution}.png',dpi=200)\n"
   ]
  },
  {
   "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(imm_clust),'leidencelltype5'] = 'immune'\n",
    "df_p.loc[df_p.leiden.isin(endo_clust),'leidencelltype5'] = 'endothelial'\n",
    "df_p.loc[df_p.leiden.isin(fb_clust),'leidencelltype5'] = 'fibroblast'\n",
    "adata.obs['leidencelltype5'] = df_p.leidencelltype5\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "figname=f'leiden__{resolution}leidencelltype5.png'\n",
    "sc.pl.umap(adata, color='leidencelltype5',ax=ax,title=figname.split('.png')[0].replace('_',' '),save=figname,wspace=0.25) #\n",
    "#fig.savefig(f'{codedir}/{s_date}/{s_sample}_umap_leiden_celltype_{resolution}.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(1.4,3.7),dpi=200)\n",
    "figname = f\"{s_sample}_{n_neighbors}neighbors_{n_markers}markers_{s_type}_leiden{resolution}.png\"\n",
    "df_p.groupby('leiden').count().loc[categories_order[::-1]].iloc[:,1].plot(kind='barh',title='Cell Count',ax=ax,width=0.7)\n",
    "fig.savefig(f'./{s_date}/barplot_{figname}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_p['slide'] = [item.split('_')[0] for item in df_p.index]\n",
    "df_p[df_p.leiden=='15'].Pan_Keratin.median()\n",
    "#df_p[df_p.leiden=='15'].Beta_catenin_c.median()\n",
    "np.quantile(df_p[df_p.leiden=='2b'].CD45,.2)\n",
    "np.quantile(df_p[df_p.leiden=='1'].Vimentin,.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### gate celltypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[df_p.CD31>.1].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='CD31')\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[df_p.CD45>2].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='CD45')\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[df_p.Beta_catenin_c>.2].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='Bcat')\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[df_p.Pan_Keratin>.5].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='panCK')\n",
    "adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[(df_p.Vimentin>1.5)].index)).replace({True:1,False:0}).values\n",
    "sc.pl.umap(adata, color='test',title='Vim')\n",
    "#adata.obs['test'] = pd.Series(adata.obs.index.isin(df_p[(df_p.FN>1)].index)).replace({True:1,False:0}).values\n",
    "#sc.pl.umap(adata, color='test',title='FN')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#manual gating\n",
    "df_p['gatedcelltype5'] = np.nan\n",
    "df_p.loc[((df_p.Vimentin>1.5)),'gatedcelltype5'] = 'fibroblast'\n",
    "df_p.loc[((df_p.Beta_catenin_c>.2) | (df_p.Pan_Keratin>.5)),'gatedcelltype5'] = 'epithelial'\n",
    "df_p.loc[df_p.CD45>2,'gatedcelltype5'] = 'immune'\n",
    "df_p.loc[df_p.CD31>.1,'gatedcelltype5'] = 'endothelial'\n",
    "df_p.gatedcelltype5 = df_p.gatedcelltype5.fillna('stromal')\n",
    "adata.obs['gatedcelltype5'] = df_p.gatedcelltype5\n",
    "fig,ax = plt.subplots(figsize=(3,2),dpi=200)\n",
    "figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}_{s_type}_gatingcelltype5.png'\n",
    "sc.pl.umap(adata, color='gatedcelltype5',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' '))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "labels = sorted(df_p.leidencelltype5.fillna('stromal').unique())\n",
    "confusion_matrix = sklearn.metrics.confusion_matrix(df_p.gatedcelltype5,df_p.leidencelltype5.fillna('stromal'),labels=labels)\n",
    "#accuracy\n",
    "i_acc = sklearn.metrics.accuracy_score(df_p.gatedcelltype5,df_p.leidencelltype5.fillna('stromal'))\n",
    "disp = sklearn.metrics.ConfusionMatrixDisplay(confusion_matrix,display_labels=labels)\n",
    "fig, ax = plt.subplots(figsize = (5,5))\n",
    "disp.plot(xticks_rotation='vertical',ax=ax) \n",
    "ax.set_ylabel('Gating Celltype')\n",
    "ax.set_xlabel(' Clustering Celltype')\n",
    "ax.set_title(f'Gating vs Clustering \\n accuracy = {i_acc:.2}')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/{s_date}/Confusion_matrix_leiden{resolution}_5cells.png',dpi=300)\n",
    "#4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.groupby('gatedcelltype5').count().CD31/len(df_p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.groupby('leidencelltype5').count().CD31/len(df_p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p['slide'] = [item.split('_')[0] for item in df_p.index]\n",
    "df_p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_out = f'{s_sample}_LeidenClusteringGating_neighbors{n_neighbors}_resolution{resolution}_markers{n_markers}_{s_type}.csv'\n",
    "df_out = df_p.merge(df_mi.loc[:,['DAPI_X','DAPI_Y']],left_index=True,right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(s_out):\n",
    "    print('saving csv')\n",
    "    df_out.to_csv(s_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize  <a name=\"viz\"></a>\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mplex_image import analyze\n",
    "df_p = pd.read_csv(f'{s_sample}_LeidenClusteringGating_neighbors{n_neighbors}_resolution{resolution}_markers{n_markers}_{s_type}.csv',index_col=0)\n",
    "df_p['slide_scene'] = [item.split('_cell')[0] for item in df_p.index]\n",
    "df_p['slide'] = [item.split('_')[0] for item in df_p.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_clust = '8'#11\n",
    "df_p[df_p.leiden==s_clust].groupby('slide').count().sort_values(by='CD31',ascending=False)[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "all_colors = mpl.cm.Set1.colors + mpl.cm.Set2.colors+ mpl.cm.Set3.colors\n",
    "#plot all groups spatially - leiden \n",
    "df_pos = analyze.celltype_to_bool(df_p,'leiden')\n",
    "df_pos.columns = [str(item) for item in df_pos.columns]\n",
    "for s_slide in sorted(set(df_p.slide)):\n",
    "    s_slide= '21'\n",
    "    fig,ax = plt.subplots(figsize=(5,3.5),dpi=200)\n",
    "    #plot negative cells\n",
    "    df_scene = df_p[(df_p.slide==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(df_pos.columns):\n",
    "        s_color=s_clust\n",
    "        s_color = str(s_color)\n",
    "        if len(df_p[(df_p.slide==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ax.scatter(data=df_p[(df_p.slide==s_slide) & (df_pos.loc[:,s_color])],x='DAPI_X',y='DAPI_Y',\n",
    "                                                                        label=f'{s_color}',s=0.1,color=all_colors[idxs])\n",
    "        break\n",
    "    ax.set_title(f\"{s_slide}\", fontsize=16)\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,ncol=2,bbox_to_anchor=(.8,.95)) \n",
    "    #fig.savefig(f'{codedir}/figures/{s_slide}_leiden{resolution}_scatterplot.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#plot all groups spatially  - manual\n",
    "df_pos = analyze.celltype_to_bool(df_p,'leiden_celltype')\n",
    "df_pos.columns = [str(item) for item in df_pos.columns]\n",
    "for s_slide in sorted(set(df_p.slide)):\n",
    "    fig,ax = plt.subplots(figsize=(5,3.5),dpi=200)\n",
    "    s_slide= '26'#'6'#'3'\n",
    "    #plot negative cells\n",
    "    df_scene =df_p[(df_p.slide==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(df_pos.columns):\n",
    "        #print(s_color)\n",
    "        #s_color = str(s_color_int)\n",
    "        if len(df_p[(df_p.slide==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ax.scatter(data=df_p[(df_p.slide==s_slide) & (df_pos.loc[:,s_color])],x='DAPI_X',y='DAPI_Y',\n",
    "                                                                        label=f'{s_color}',s=0.1,color=mpl.cm.tab10.colors[idxs])\n",
    "        #break\n",
    "    ax.set_title(f\"{s_slide}\", fontsize=16)\n",
    "    ax.axis('equal')\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "    for spines in ['left','right','top','bottom']:\n",
    "        ax.spines[spines].set_visible(False)\n",
    "\n",
    "    #break\n",
    "    plt.legend(markerscale=10,framealpha=.1,bbox_to_anchor=(.85,.5)) \n",
    "    fig.savefig(f'{codedir}/figures/{s_slide}_gated_cells_scatterplot.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#old clusters (based on the )\n",
    "'''\n",
    "d_group = {'10': 'stroma_rich', '11': 'stroma_rich', '12': 'stroma_rich', '13': 'stroma_rich', '14': 'stroma_poor',\n",
    " '15': 'stroma_poor', '16': 'stroma_rich', '17': 'stroma_rich', '18': 'stroma_poor', '19': 'stroma_rich', '2': 'stroma_rich',\n",
    " '20': 'stroma_rich', '21': 'stroma_poor', '23': 'stroma_poor', '24': 'stroma_poor', '25': 'stroma_rich',\n",
    " '26': 'stroma_poor', '27': 'stroma_poor', '28': 'stroma_rich', '29': 'stroma_rich', '3': 'stroma_rich',\n",
    " '31': 'stroma_poor', '32': 'stroma_rich', '33': 'stroma_rich', '34': 'stroma_rich', '35': 'stroma_rich', '36': 'stroma_rich',\n",
    " '37': 'stroma_rich', '39': 'stroma_rich', '4': 'stroma_rich', '40': 'stroma_rich', '41': 'stroma_rich', '5': 'stroma_rich',\n",
    " '6': 'stroma_rich', '7': 'stroma_poor', '8': 'stroma_rich', '9': 'stroma_rich'}\n",
    "d_group3 = {'10': 'imm', '11': 'str', '12': 'imm', '13': 'imm', '14': 'epi', '15': 'epi', '16': 'imm', '17': 'imm', '18': 'epi',\n",
    " '19': 'str', '2': 'imm', '20': 'str', '21': 'epi', '23': 'epi', '24': 'epi', '25': 'str', '26': 'epi', '27': 'epi', '28': 'imm',\n",
    " '29': 'imm', '3': 'imm', '31': 'epi', '32': 'imm', '33': 'str', '34': 'str', '35': 'imm', '36': 'imm', '37': 'imm', '39': 'imm',\n",
    " '4': 'imm', '40': 'imm', '41': 'imm', '5': 'imm', '6': 'str', '7': 'epi', '8': 'str', '9': 'imm'}\n",
    "\n",
    "#new\n",
    "d_group = {'1': 'stroma_rich', '10': 'stroma_rich', '11': 'stroma_rich', '12': 'stroma_rich', '13': 'stroma_rich',\n",
    " '14': 'epithelial_rich', '15': 'epithelial_rich', '16': 'stroma_rich', '17': 'stroma_rich', '18': 'epithelial_rich',\n",
    " '19': 'epithelial_rich', '2': 'epithelial_rich', '20': 'stroma_rich', '21': 'epithelial_rich',\n",
    " '22': 'epithelial_rich', '23': 'epithelial_rich', '24': 'epithelial_rich', '25': 'epithelial_rich',\n",
    " '26': 'epithelial_rich', '27': 'epithelial_rich', '28': 'stroma_rich', '29': 'stroma_rich', '3': 'stroma_rich',\n",
    " '30': 'stroma_rich', '31': 'epithelial_rich', '32': 'stroma_rich', '33': 'epithelial_rich', '34': 'stroma_rich',\n",
    " '35': 'stroma_rich', '36': 'stroma_rich', '37': 'stroma_rich', '38': 'epithelial_rich', '39': 'stroma_rich',\n",
    " '4': 'stroma_rich', '40': 'stroma_rich', '41': 'stroma_rich', '5': 'stroma_rich', '6': 'stroma_rich',\n",
    " '7': 'epithelial_rich', '8': 'epithelial_rich', '9': 'stroma_rich'}\n",
    "#\n",
    "\n",
    "df_p['Group'] = df_p.slide.map(d_group)\n",
    "se_plot = df_p[df_p.Group!='not_tumor'].groupby(['leiden','Group']).count().loc[categories_order[::-1]].iloc[:,1]\n",
    "se_plot_frac = se_plot/df_p[df_p.Group!='not_tumor'].groupby(['Group']).count().iloc[:,1]\n",
    "fig,ax=plt.subplots(figsize=(2.5,5),dpi=200)\n",
    "sns.barplot(data=se_plot_frac.reset_index(),y='leiden',x='Beta_catenin_n',hue='Group',orient='h',order=categories_order,ax=ax,hue_order = ['stroma_rich','epithelial_rich'])\n",
    "ax.set_xlabel('Fraction of Cells')\n",
    "ax.legend(bbox_to_anchor=(0.3,1.01))\n",
    "fig.savefig(f'./{s_date}/barplot_{n_neighbors}neighbors_{len(marker_genes)}markers_leiden{resolution}_bygroup.png')\n",
    "#\n",
    "df_p['slide'] = [item.split('_')[0] for item in df_p.index]\n",
    "df_p['Group'] = df_p.slide.map(d_group3)\n",
    "se_plot = df_p[df_p.Group!='not_tumor'].groupby(['leiden','Group']).count().loc[categories_order[::-1]].iloc[:,1]\n",
    "se_plot_frac = se_plot/df_p[df_p.Group!='not_tumor'].groupby(['Group']).count().iloc[:,1]\n",
    "fig,ax=plt.subplots(figsize=(2.5,5),dpi=200)\n",
    "sns.barplot(data=se_plot_frac.reset_index(),y='leiden',x='Beta_catenin_n',hue='Group',orient='h',order=categories_order,ax=ax)\n",
    "ax.set_xlabel('Fraction of Cells')\n",
    "fig.savefig(f'./{s_date}/barplot_{n_neighbors}neighbors_{len(marker_genes)}markers_leiden{resolution}_bygroup3.png')\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# subtyping <a name=\"subtype\"></a>\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '20211109_MIBI'\n",
    "resolution = 0.6\n",
    "i_cells = '' #'5' # \n",
    "df_p = pd.read_csv(f'results_{s_sample}_leiden{resolution}_gated_celltype{i_cells}.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p['slide'] = [item.split('_')[0] for item in df_p.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "s_type = 'gated_celltype' #'leiden' #\n",
    "df = df_p.loc[:,['CD31',s_type,'slide']].groupby([s_type,'slide']).CD31.count()/df_p.loc[:,['CD31','slide']].groupby(['slide']).CD31.count()\n",
    "df_prop = df.unstack().T.fillna(0)\n",
    "ls_col = df_prop.columns.tolist()\n",
    "adata = sc.AnnData(df_prop)\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)}_{s_type}.png')\n",
    "plt.tight_layout()\n",
    "sc.pl.pca_variance_ratio(adata,n_pcs=32, log=True)\n",
    "\n",
    "#umap\n",
    "for n_neighbors in [5]: #,30\n",
    "    results_file = f'{s_date}/{s_sample}_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.h5ad'\n",
    "    print('calc umap')\n",
    "    # calculate neighbors \n",
    "    sc.pp.neighbors(adata, n_neighbors=n_neighbors) #, method='rapids'\n",
    "    sc.tl.umap(adata)\n",
    "\n",
    "    #color by markers\n",
    "    \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.5',ncols=4,save=figname,size=250)\n",
    "    #survival\n",
    "    '''\n",
    "    adata.obs['slide'] = adata.obs.index.tolist()\n",
    "    adata.obs['Group'] = adata.obs.slide.map(d_a).fillna('unknown')\n",
    "    #adata.obs['Group'] = adata.obs.index.map(d_a).fillna('unknown')\n",
    "            \n",
    "    figname = f\"Umap_Group_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.png\"\n",
    "    title=figname.split('.png')[0].replace('_',' ')\n",
    "    fig,ax = plt.subplots(figsize=(3,2), dpi=200)\n",
    "    sc.pl.umap(adata, color='Group',wspace=.25,save=figname,ax=ax,title=title,size=50)\n",
    "    \n",
    "    # time\n",
    "    figname = f\"Umap_survival_{n_neighbors}neighbors_{len(ls_col)}markers_{s_type}.png\"\n",
    "    title=figname.split('.png')[0].replace('_',' ')\n",
    "    fig,ax = plt.subplots(figsize=(3,2), dpi=200)\n",
    "    adata.obs['survival'] =adata.obs.slide.map(dict(zip(df_a.index.astype('str'),df_a.Survival_time)))\n",
    "    sc.pl.umap(adata, color='survival',size=50,ax=ax,title=title,save=figname)\n",
    "    '''\n",
    "    #save results\n",
    "    if not os.path.exists(results_file):\n",
    "        adata.write(results_file)\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#leiden\n",
    "for resolution in [0.2]: #0.3,\n",
    "    results_file = f'{s_sample}_{n_neighbors}neighbors_{n_markers}markers_{s_type}_leiden{resolution}.h5ad'\n",
    "    sc.tl.leiden(adata,resolution=resolution)\n",
    "    adata.write(results_file)\n",
    "    fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "    figname=f'leiden_{n_neighbors}_{n_markers}_{resolution}.png'\n",
    "    sc.pl.umap(adata, color='leiden',ax=ax,title=figname.split('.png')[0].replace('_',' '),wspace=.25,save=figname,size=40)\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#sns.clustermap(df_prop)\n",
    "\n",
    "d_a = dict(zip(adata.obs.index,adata.obs.leiden))\n",
    "d_color = dict(zip(adata.obs.leiden.unique(),sns.color_palette()[0:len(adata.obs.leiden.unique())]))\n",
    "row_colors = df_prop.index.map(d_a).map(d_color)\n",
    "g = sns.clustermap(df_prop,figsize=(5,15),cmap='viridis',\n",
    "        row_colors=row_colors)\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.savefig(f'{codedir}/{s_date}/{s_sample}_clustermap_Group.png',dpi=200)\n",
    "marker_genes = df_prop.iloc[:,g.dendrogram_col.reordered_ind].columns.tolist()\n",
    "categories_order = df_prop.iloc[g.dendrogram_row.reordered_ind,:].index.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if resolution ==0.2:\n",
    "    ls_epi = ['15', '22', '18', '24', '7', '14', '21', '23', '27', '26', '31',]#1\n",
    "    ls_str = ['38', '11', '20', '34', '6', '25', '33', '19', '8',] #2\n",
    "    ls_imm = [ '32', '35', '16', '28', '2', '36', 't1', '3', '5', '12', '9', '17', '39', '29', '37', '41', '10', '40', '13', '4'] #0\n",
    "    ls_myl = []\n",
    "elif resolution == 0.7:\n",
    "    ls_epi = ['15', '22', '18', '24', '7', '14', '21', '23', '27', '26', '31',]\n",
    "\n",
    "    ls_str = [ '38', '11', '20', '34',\n",
    "     '6', '25', '33', '19', '8','32', '35', '16', '28', '2', '36', '1', '3', '5', '12', '9', '17', '39', '29', '37',\n",
    "     '41', '10', '40', '13', '4']\n",
    "    ls_imm = ['10', '12', '13', '16', '17', '2', '28', '29', '3', '32', '35', '36',\n",
    "           '37', '39', '4', '40', '41', '5', '9', 't1']\n",
    "    ls_myl = []\n",
    "    ls_lym = []\n",
    "elif i_Cells==5:\n",
    "    ls_epi =  ['15', '18', '22', '24','14', '21', '23', '26', '27', '31', '7']\n",
    "    ls_myl = ['12', '17', '3', '39', '5', '9', 't1']\n",
    "    ls_lym = ['16', '28', '32', '35']\n",
    "    ls_str = ['10', '13', '2', '29', '36', '37', '4', '40', '41','11', '19', '20', '25', '33', '34', '38', '6', '8',]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clustermap\n",
    "df_p = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "df_p['leiden'] = adata.obs['leiden']\n",
    "d_replace = {}\n",
    "df_plot = df_p.groupby('leiden').mean()\n",
    "df_plot.index.name = f'leiden {resolution}'\n",
    "g = sns.clustermap(df_plot,z_score=1,figsize=(4,3.5),cmap='viridis',\n",
    "                   vmin=-1.5,vmax=1.5)\n",
    "#g.fig.suptitle(f'leiden {resolution}',x=.9) \n",
    "g.savefig(f'{codedir}/{s_date}/{s_sample}_clustermap_{s_type}_subtype_{resolution}.png',dpi=200)\n",
    "marker_genes = df_p.groupby('leiden').mean().iloc[:,g.dendrogram_col.reordered_ind].columns.tolist()\n",
    "categories_order = df_p.groupby('leiden').mean().iloc[g.dendrogram_row.reordered_ind,:].index.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(2.5,3),dpi=200)\n",
    "df_p.groupby('leiden').count().loc[categories_order[::-1]].iloc[:,1].plot(kind='barh',title='Patient Count',ax=ax)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./{s_date}/barplot_{s_sample}_clustermap_leiden_subtype_{resolution}.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p[df_p.leiden=='0'].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "categories_order\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# survival <a name=\"surv\"></a>\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lifelines\n",
    "from lifelines import KaplanMeierFitter\n",
    "from lifelines.statistics import multivariate_logrank_test\n",
    "from lifelines import CoxPHFitter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "df_a.index = df_a.index.astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df_p.merge(df_a.loc[:,['Survival','Survival_time']],left_index=True,right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = df['Survival_time']     ## time to event\n",
    "E = df['Survival']      ## event occurred or censored\n",
    "groups = df.loc[:,'leiden']  \n",
    "kmf1 = KaplanMeierFitter() ## instantiate the class to create an object\n",
    "fig, ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "for s_group in sorted(groups.unique()):\n",
    "    i1 = (groups == s_group)\n",
    "    kmf1.fit(T[i1], E[i1], label=s_group)    ## fit thedata\n",
    "    kmf1.plot(ax=ax,ci_show=False)\n",
    "    #print(kmf1.median_survival_time_)\n",
    "results = multivariate_logrank_test(event_durations=T, groups=groups, event_observed=E)\n",
    "ax.set_title(f'{s_type} {resolution} \\n p={results.summary.p[0]:.2}')\n",
    "ax.legend(loc='upper right')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{s_date}/KM_{resolution}_{s_type}.png',dpi=200)\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[df.index.isin(ls_epi),'Group'] ='stroma_poor'\n",
    "df.loc[df.index.isin(ls_str),'Group'] ='stroma_rich'\n",
    "if len(ls_imm) > 0:\n",
    "    df.loc[df.index.isin(ls_imm),'Group'] ='stroma_rich'    \n",
    "if len(ls_myl) > 0:\n",
    "    df.loc[df.index.isin(ls_myl),'Group'] ='myl'\n",
    "    df.loc[df.index.isin(ls_lym),'Group'] ='lym'\n",
    "    \n",
    "    \n",
    "\n",
    "s_type = 'Group'\n",
    "groups = df.loc[:,'Group']  \n",
    "kmf1 = KaplanMeierFitter() ## instantiate the class to create an object\n",
    "fig, ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "for s_group in sorted(groups.unique()):\n",
    "    i1 = (groups == s_group)\n",
    "    kmf1.fit(T[i1], E[i1], label=s_group)    ## fit thedata\n",
    "    kmf1.plot(ax=ax,ci_show=False)\n",
    "    #print(kmf1.median_survival_time_)\n",
    "results = multivariate_logrank_test(event_durations=T, groups=groups, event_observed=E)\n",
    "ax.set_title(f'{s_type} {resolution} \\n p={results.summary.p[0]:.2}')\n",
    "ax.legend(loc='upper right')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{s_date}/KM_{resolution}_{s_type}.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dict(zip(df.index,df.Group3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[df.index.isin(ls_epi),'Group3'] ='epi'\n",
    "df.loc[df.index.isin(ls_str),'Group3'] ='str'\n",
    "df.loc[df.index.isin(ls_imm),'Group3'] ='imm'\n",
    "if len(ls_myl) > 0:\n",
    "    df.loc[df.index.isin(ls_myl),'Group3'] ='myl'\n",
    "    df.loc[df.index.isin(ls_lym),'Group3'] ='lym'\n",
    "    \n",
    "    \n",
    "\n",
    "s_type = 'Group3'\n",
    "groups = df.loc[:,'Group3']  \n",
    "kmf1 = KaplanMeierFitter() ## instantiate the class to create an object\n",
    "fig, ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "for s_group in sorted(groups.unique()):\n",
    "    i1 = (groups == s_group)\n",
    "    kmf1.fit(T[i1], E[i1], label=s_group)    ## fit thedata\n",
    "    kmf1.plot(ax=ax,ci_show=False)\n",
    "    #print(kmf1.median_survival_time_)\n",
    "results = multivariate_logrank_test(event_durations=T, groups=groups, event_observed=E)\n",
    "ax.set_title(f'{s_type} {resolution} \\n p={results.summary.p[0]:.2}')\n",
    "ax.legend(loc='upper right')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{s_date}/KM_{resolution}_{s_type}.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cph = CoxPHFitter()   ## Instantiate the class to create a cph object\n",
    "cph.fit(df.drop(['leiden','Group','immune','stromal','endothelial'],axis=1), 'Survival_time', event_col='Survival')   ## Fit the data to train the model\n",
    "#cph.print_summary()    ## HAve a look at the significance of the features\n",
    "fig, ax = plt.subplots()\n",
    "cph.plot(ax=ax)\n",
    "ax.set_title('Cox Proportional Hazards')\n",
    "plt.tight_layout()\n",
    "fig.savefig('Cox_PH_epithelial.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "df['es_ratio'] = df.loc[:,'epithelial']/df.loc[:,['endothelial','immune','stromal']].sum(axis=1)\n",
    "cph = CoxPHFitter()  \n",
    "cph.fit(df.drop(['leiden','Group','immune','stromal','endothelial','epithelial'],axis=1), 'Survival_time', event_col='Survival')   ## Fit the data to train the model\n",
    "#cph.print_summary()    ## HAve a look at the significance of the features\n",
    "fig, ax = plt.subplots()\n",
    "cph.plot(ax=ax)\n",
    "ax.set_title('Cox Proportional Hazards')\n",
    "plt.tight_layout()\n",
    "fig.savefig('Cox_PH_es_ratio.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#single var"
   ]
  },
  {
   "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
}