20211109_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/rasp-mibi'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change to correct directory\n",
    "datadir =  f'{codedir}/rawdata'\n",
    "s_date = '20211206'#'20211109'#'20211029'\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",
    "1. [Load manual](#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 Gated Celltypes <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",
    "df_pos = pd.DataFrame()\n",
    "for s_file in os.listdir('intermediate_data/protein_positivity'):\n",
    "    if s_file.find('.csv') > -1:\n",
    "        df=pd.read_csv(f'./intermediate_data/protein_positivity/{s_file}')\n",
    "        df['Sample'] = s_file.split('.')[0]\n",
    "        df_pos = df_pos.append(df)\n",
    "    #break\n",
    "df_pos.index = df_pos.Sample.astype('str') + '_' + df_pos.index.astype('str')\n",
    "df_pos = df_pos.rename(d_name,axis=1)"
   ]
  },
  {
   "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": [
    "#mean intensity and centroid XY\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",
    "d_name = dict(zip(df_name.index.astype('str'),df_name.Biomarker))\n",
    "\n",
    "df_mi_mean = pd.read_csv(f'./intermediate_data/proteinexpression_df.csv',index_col=0)\n",
    "\n",
    "df_mi = pd.DataFrame()\n",
    "for s_file in os.listdir('intermediate_data/protein_expression'):\n",
    "    df=pd.read_csv(f'./intermediate_data/protein_expression/{s_file}')\n",
    "    df['Sample'] = s_file.split('.')[0]\n",
    "    df_mi = df_mi.append(df)\n",
    "    #break\n",
    "df_mi.index = df_mi.Sample.astype('str') + '_' + df_mi.index.astype('str')\n",
    "df_mi = df_mi.rename(d_name,axis=1)\n",
    "df_xy = df_mi.loc[:,['Row','Column','Sample']].rename({'Row':'DAPI_Y','Column':'DAPI_X'},axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_name[df_name.Purpose=='Lineage'].Biomarker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_name.Biomarker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_shared = ['CD11b', 'CD11c','CD31','CD45','CD8','dsDNA','Beta catenin','FoxP3','Keratin6','Keratin17', 'Ki67', 'Pan-Keratin','Vimentin','SMA']"
   ]
  },
  {
   "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",
    "#df_mi.columns = [item.split('_'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_drop = ['Au', 'Background', 'Ca', 'Celltype', 'Column', 'Fe', 'Na', 'P','Row','Sample' ,'Si', 'Ta']"
   ]
  },
  {
   "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_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_mi.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_exclude = pd.DataFrame()\n",
    "for s_file in os.listdir():\n",
    "    if s_file.find('exclude_RS') > -1:\n",
    "        print(s_file)\n",
    "        df_exclude = df_exclude.append(pd.read_csv(s_file,index_col=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_filter_mi = df_mi.merge(df_xy.loc[:,['DAPI_X', 'DAPI_Y']], left_index=True,right_index=True)\n",
    "#df_filter_mi.rename({'slidescene':'slide_scene'},axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_filter_mi#.head()\n",
    "#s_date\n",
    "s_sample='MIBI'"
   ]
  },
  {
   "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": true
   },
   "outputs": [],
   "source": [
    "s_sample = '20211109_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', 'CD11b', 'CD11c', 'CD138', 'CD16', 'CD20', 'CD209',\n",
    "       'CD3', 'CD31', 'CD4', 'CD45', 'CD45RO', 'CD56', 'CD63', 'CD68', 'CD8',\n",
    "       'dsDNA', 'EGFR', 'FoxP3', 'HLA-DR',  'IDO', # 'H3K27me3', 'H3K9ac', 'phospho-S6','HLA',\n",
    "       'Keratin17', 'Keratin6', 'Ki67', 'Lag3', 'MPO', 'p53', 'Pan-Keratin',\n",
    "       'PD-L1', 'PD1', 'SMA', 'Vimentin'],\n",
    "    ['CD11b', 'CD11c','CD31','CD45','CD8','dsDNA',\n",
    "       'Beta catenin','FoxP3','Keratin6','Keratin17', 'Ki67', 'Pan-Keratin','Vimentin','SMA'], #not good clustering\n",
    "\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",
    "    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",
    "\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 [15]: #,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.replace({'t1':'1'}).unique()"
   ]
  },
  {
   "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 = '20211109_MIBI'\n",
    "s_type = 'Filtered'\n",
    "s_date = s_sample.split('_')[0]\n",
    "n_neighbors = 15\n",
    "n_markers = 32\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.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 = '20211109_MIBI'\n",
    "s_type = 'Filtered'\n",
    "s_date = s_sample.split('_')[0]\n",
    "n_neighbors = 16\n",
    "n_markers = 32\n",
    "resolution = 0.6\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": {},
   "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']\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 resolution == 0.5:\n",
    "    tum_clust = [ '15', '13', '16', '5', '8','9','2','6','14','7', '20',]\n",
    "    str_clust = ['0', '1', '10', '11', '12', '17', '18', '19', '21', '3', '4']\n",
    "    imm_clust = [] \n",
    "    endo_clust = []\n",
    "if resolution == 0.6:\n",
    "    imm_clust = ['7', '1', '23', '14', '0', '3']\n",
    "    endo_clust = ['13','23']#[]\n",
    "    myl_clust = ['20']\n",
    "    tum_clust = ['19', '16', '22', '24',  '12', '8', '18','4', '10', '11', '6', '2', '9'] #'17', '21',\n",
    "    str_clust = ['17', '21','15', '5','7', '23',  '14', '1', '0', '3', '20','13']\n",
    "df_p.groupby('leiden').mean().CD11b.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": [
    "#clustermap\n",
    "#df_p = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "d_replace = {}\n",
    "d_replace.update(dict(zip(tum_clust,[f'{item}\\nep' for item in tum_clust])))\n",
    "d_replace.update(dict(zip(str_clust,[f'{item}\\ns' for item in str_clust])))\n",
    "d_replace.update(dict(zip(imm_clust,[f'{item}\\nim' for item in imm_clust])))\n",
    "d_replace.update(dict(zip(endo_clust,[f'{item}\\nen' for item in endo_clust])))\n",
    "d_replace.update(dict(zip(myl_clust,[f'{item}\\nmyl' for item in myl_clust])))\n",
    "df_plot = df_p.groupby('leiden').mean().rename(d_replace)\n",
    "df_plot.index.name = f'leiden {resolution}'\n",
    "g = sns.clustermap(df_plot,z_score=1,figsize=(10,10),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_leiden_{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": [
    "#categories_order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(2.5,5),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)\n",
    "fig.savefig(f'./{s_date}/barplot_{figname}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "df_p['slide'] = [item.split('_')[0] for item in df_p.index]\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='CD11b',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}_bygroup.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\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='CD11b',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')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.loc[df_p.leiden.isin(tum_clust),'leiden_celltype'] = 'epithelial'\n",
    "df_p.loc[df_p.leiden.isin(str_clust),'leiden_celltype'] = 'stromal'\n",
    "#df_p.loc[df_p.leiden.isin(imm_clust),'leiden_celltype'] = 'immune'\n",
    "#df_p.loc[df_p.leiden.isin(endo_clust),'leiden_celltype'] = 'endothelial'\n",
    "#df_p.loc[df_p.leiden.isin(myl_clust),'leiden_celltype'] = 'myeloid'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adata.obs['leiden_celltype'] = df_p.leiden_celltype\n",
    "fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "figname=f'leiden__{resolution}_leiden_celltype.png'\n",
    "sc.pl.umap(adata, color='leiden_celltype',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": "markdown",
   "metadata": {},
   "source": [
    "# gate stromal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p['gated_celltype'] = df_p.leiden_celltype\n",
    "i_thresh_cd45 = 0.75\n",
    "i_thresh_cd31 = 0.25\n",
    "i_thresh_cd11b = 0.25\n",
    "#gated cells\n",
    "df_p.loc[((df_p.leiden_celltype == 'stromal') & (df_p.CD31 > i_thresh_cd31)),'gated_celltype'] = 'endothelial'\n",
    "df_p.loc[((df_p.leiden_celltype == 'stromal') & (df_p.CD45 > i_thresh_cd45) & (df_p.CD31 <= i_thresh_cd31)),'gated_celltype'] = 'immune'\n",
    "df_p.loc[((df_p.gated_celltype=='immune') & (df_p.CD11b >= i_thresh_cd11b)),'gated_celltype'] = 'myeloid'\n",
    "df_p.loc[((df_p.gated_celltype=='immune') & (df_p.CD11b < i_thresh_cd11b)),'gated_celltype'] = 'lymphoid'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata.obs['gated_celltype'] = df_p.gated_celltype\n",
    "fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "figname=f'_gated_celltype.png'\n",
    "sc.pl.umap(adata, color='gated_celltype',ax=ax,save=figname,title=figname.split('.png')[0].replace('_',' ')) #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi.rename({'Row':'DAPI_Y','Column':'DAPI_X'},axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.merge(df_mi.rename({'Row':'DAPI_Y','Column':'DAPI_X'},axis=1).loc[:,['DAPI_X','DAPI_Y']],left_index=True,right_index=True).to_csv(f'results_{s_sample}_leiden{resolution}_gated_celltype5.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sub Cluster  <a name=\"kmeans\"></a>\n",
    "\n",
    "cluster patients by celltype\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize  <a name=\"viz\"></a>\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '20211109_MIBI'\n",
    "resolution = 0.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mplex_image import analyze\n",
    "df_p = pd.read_csv(f'results_{s_sample}_leiden{resolution}_gated_celltype.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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#plot all groups spatially  - manual\n",
    "df_pos = analyze.celltype_to_bool(df_p,'gated_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= '3'#'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.set_xticklabels('')\n",
    "    #ax.set_yticklabels('')\n",
    "    #break\n",
    "    plt.legend(markerscale=10,framealpha=.1) \n",
    "    fig.savefig(f'{codedir}/figures/{s_slide}_gated_cells_scatterplot.png')\n",
    "    break"
   ]
  },
  {
   "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= '6'\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 = 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": [
    "#sorted(df_p.slide.unique()) missing 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# subtyping <a name=\"subtype\"></a>\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "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
}