20211027_RS-TMA_cluster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster Mouse Samples\n",
    "\n",
    "**Samples:** \n",
    "- RS-TMA-5\n",
    "\n",
    "**Method**: Cluster Umap with Leiden algorithm \n",
    "\n",
    "Compare to images\n",
    "\n",
    "compare to histological 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/20210000/cmIF_2021-07-07_RS-mTMA'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change to correct directory\n",
    "datadir =  f'{codedir}'\n",
    "s_date = '20211110'#'20211029'\n",
    "os.chdir(codedir)\n",
    "if not os.path.exists(s_date):\n",
    "    os.mkdir(s_date)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of contents <a name=\"contents\"></a>\n",
    "1. [Load 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)"
   ]
  },
  {
   "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": [
    "s_out = '20211109_RS-mTMA-5'\n",
    "df_man = pd.read_csv(f'{s_out}_GatedPositiveCellNames.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_man.index = df_man.jindex\n",
    "#df_man.index = [item.replace('ad60','2021.02.03-adeno60').replace('ad82','2021.02.04-adeno82') for item in df_man.index]"
   ]
  },
  {
   "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": [
    "#sample\n",
    "ls_sample = ['RS-mTMA-5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mean intensity and centroid XY\n",
    "df_xy = pd.DataFrame()\n",
    "df_mi = pd.DataFrame()\n",
    "for s_sample in ls_sample:\n",
    "    #1027\n",
    "    #df_mi = df_mi.append(pd.read_csv(f'{datadir}/features_{s_sample}_patched_MeanIntensity_Shape_DAPI1_DAPI8_subtractedregisteredimages.csv',index_col=0))\n",
    "    #1029\n",
    "    df_mi = df_mi.append(pd.read_csv(f'{datadir}/features_{s_sample}_patched_MeanIntensity_Shape__subtractedregisteredimages.csv',index_col=0))\n",
    "    df_xy = df_xy.append(pd.read_csv(f'{datadir}/features_{s_sample}_CentroidXY.csv',index_col=0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sorted(df_mi.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_mi.columns = [item.split('_')[0] for item in df_mi.columns]\n",
    "ls_drop = ['DAPI0_nuclei', 'DAPI1_nuclei', 'DAPI2_nuclei', 'DAPI3_nuclei', 'DAPI4_nuclei', 'DAPI5Q_nuclei', 'DAPI5_nuclei',\n",
    " 'DAPI6_nuclei', 'DAPI7_nuclei', #'DAPI8_nuclei','R1c2_nuclei',(AF)\n",
    "           'R0c2_nuclei', 'R0c2_perinuc5', 'R0c3_nuclei', 'R0c3_perinuc5', 'R0c4_nuclei',\n",
    " 'R0c4_perinuc5', 'R0c5_nuclei', 'R0c5_perinuc5',  'R1c2_perinuc5', 'R4c2_nuclei', 'R4c2_perinuc5', 'R5Qc2_nuclei',\n",
    " 'R5Qc2_perinuc5', 'R5Qc3_nuclei', 'R5Qc3_perinuc5', 'R5Qc4_nuclei', 'R5Qc4_perinuc5', 'R5Qc5_nuclei', 'R5Qc5_perinuc5','cell_area',\n",
    " 'cell_eccentricity', 'cell_euler', 'cytoplasm_area', 'cytoplasm_eccentricity', 'cytoplasm_euler','Ecad_negative',\n",
    " 'EpCAM_cellmem2p25','Ecad_cellmem2p25',\n",
    " #failed\n",
    " 'CC3-488_nuclei', 'CC3-647_nuclei', 'CD4_perinuc5','CD3_perinuc5','F4-80_perinuc5','MHCII_cytoplasm','ELN_perinuc5','NFkB_cytoplasm',\n",
    " 'Arg1_cytoplasm','panCK_cytoplasm',\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_mi[(~df_mi.index.isin(df_exclude.index))]\n",
    "#df_mi[(df_mi.DAPI8 > 1000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_filter_mi = df_mi[(df_mi.DAPI8 > 1000) & (~df_mi.index.isin(df_exclude.index))].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"
   ]
  },
  {
   "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": [
    "os.chdir(codedir)\n",
    "s_date = '20211027' #326562 rows\n",
    "s_date = '20211029' #748728  rows\n",
    "s_date = '20211109' #722760 rows\n",
    "s_date = '20211110' #718696 rows 16 neighbors\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": [
    "'''\n",
    "b_control = ((df_both.index.str.contains('JE-TMA')) & (df_both.slide_scene!='JE-TMA-42_scene12')) # & (df_mi.loc[:,ls_pos].isna().sum(axis=1)==0)\n",
    "#fit\n",
    "data = df_both.loc[b_control,ls_pos].T\n",
    "batch = df_both.loc[b_control,'batch']\n",
    "gamma_star, delta_star, stand_mean, var_pooled = normalize.combat_fit(data, batch)\n",
    "#transform\n",
    "#data = df_both.loc[:,ls_pos].T\n",
    "#batch = df_both.loc[:,'batch']\n",
    "#bayesdata = normalize.combat_transform(data,batch,gamma_star,delta_star,stand_mean, var_pooled)\n",
    "#df_norm = bayesdata.T\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)"
   ]
  },
  {
   "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",
    "import re\n",
    "df_a = pd.read_csv('RS-MYC-PTENfl_TMA_subgroups.csv')\n",
    "df_a.loc[~df_a.GroupA.isna(),'Groups'] = 'stroma_rich'\n",
    "df_a.loc[~df_a.GroupB.isna(),'Groups'] = 'stroma_poor'\n",
    "df_a.loc[df_a.Groups.isna(),'Groups'] = 'not_tumor'\n",
    "df_a['Xiaoyan'] = df_a.Xiaoyan.fillna(df_a.TMA)\n",
    "df_a.index = df_a.TMA #drop(['oldE1','oldF6'],inplace=True)\n",
    "d_a = dict(zip(df_a.index,df_a.Groups))\n",
    "df_a['slide_scene'] = [f'RS-mTMA-5_scene{re.sub(r\"[0-9]+\", \"\", item)}{int(re.sub(r\"[A-Z]+\", \"\", item)):02d}' for item in df_a.index]\n",
    "df_a['Group'] = df_a.Xiaoyan.map(d_a).fillna('not_tumor')      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_a.loc['A1']\n",
    "#df_a.loc['B1']\n",
    "#df_a.loc['C1']\n",
    "#df_a.loc['D1']\n",
    "#df_a.loc['E1'] 20211110\n",
    "#df_a.loc['G6']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "s_sample = '20211110_RS-mTMA-5'\n",
    "b_scale=True\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 = [[ 'Ki67', 'LamAC', 'SOX9',#'FoxP3',\n",
    "       'pMYC', 'CD103', 'CD11b', 'CD11c', 'CD31', 'CD45', 'CD8', 'ColIV',\n",
    "       'ColVI', 'CK5', 'Ecad', 'EpCAM', 'Gal3', 'S100A6', 'Vim', 'aSMA',\n",
    "       'area', 'eccentricity','DAPI8','R1c2'],\n",
    "     [ 'Ki67', 'LamAC', 'SOX9', #'FoxP3','CD103','CD11c', 'CD31','CD8',\n",
    "       'pMYC',  'CD11b',  'CD45',  'ColIV',\n",
    "       'ColVI', 'CK5', 'Ecad', 'EpCAM', 'Gal3', 'S100A6', 'Vim', 'aSMA',\n",
    "       'area', 'eccentricity']]\n",
    "\n",
    "for s_type in ['BGSubtracted','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",
    "        break\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('_cell')[0] for item in adata.obs.index]\n",
    "        if b_scale:\n",
    "            sc.pp.scale(adata, zero_center=False, max_value=20)\n",
    "            s_sample = s_sample + '_s'\n",
    "        #load annotation\n",
    "        # mapping\n",
    "        d_a = dict(zip(df_a.slide_scene,df_a.Group))\n",
    "        adata.obs['Group'] = adata.obs.slide.map(d_a)\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",
    "                # mapping\n",
    "                d_a = dict(zip(df_a.slide_scene,df_a.Group))\n",
    "                adata.obs['Group'] = adata.obs.slide.map(d_a)\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",
    "            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",
    "            break\n",
    "        break\n",
    "    break\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#screen all\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",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# weird ones (fixed)\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)"
   ]
  },
  {
   "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": [
    "len(ls_col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '20211110_RS-mTMA-5_s'\n",
    "s_type = 'BGSubtracted'\n",
    "s_date = s_sample.split('_')[0]\n",
    "n_neighbors = 15\n",
    "n_markers = 23 #24\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.4, 0.5,0.6]: #0.5,\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 ='20211110_RS-mTMA-5_s'\n",
    "s_type = 'BGSubtracted'\n",
    "s_date = s_sample.split('_')[0]\n",
    "n_neighbors = 15#16\n",
    "n_markers = 24\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": [
    "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['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",
    "b_annot = False\n",
    "if results_file=='20211110_RS-mTMA-5_s_15neighbors_24markers_BGSubtracted_leiden0.6.h5ad':\n",
    "        tum_clust = [] \n",
    "        imm_clust = []\n",
    "        str_clust = []\n",
    "        endo_clust = []\n",
    "    \n",
    "elif resolution == 0.4:  \n",
    "    if n_neighbors == 15:\n",
    "        tum_clust = ['3', '5', '10', '0', '14', '11', '6', '17',   '13'] \n",
    "        imm_clust = ['14']\n",
    "        str_clust = ['1', '16', '2', '4', '9','12','8']\n",
    "        endo_clust = ['15','7']\n",
    "elif resolution == 0.5:\n",
    "    if n_neighbors == 16:\n",
    "        tum_clust = ['7', '1',  '5', '8', '2', '11', '16','12'] #15\n",
    "        str_clust = ['0', '10',   '3','16','14','6','9','13','4','15','18'] #'12',\n",
    "        imm_clust = ['16','14','6'] # []\n",
    "        endo_clust = ['17','9'] # []\n",
    "    elif n_neighbors == 15:\n",
    "        tum_clust = ['1', '4', '0', '15', '11', '19', '9', '5', '22', '18','13','6','21']\n",
    "        str_clust = ['10', '12', '16', '17', '2', '3',  '7', '8']\n",
    "        imm_clust = ['19'] \n",
    "        endo_clust = ['14','20']\n",
    "else:\n",
    "    print('not annot')\n",
    "    b_annot = False\n",
    "\n",
    "df_p.groupby('leiden').mean().CD45.sort_values(ascending=False)#.index\n",
    "#set(df_p.groupby('leiden').mean().area.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",
    "#df_p['leiden'] = adata.obs['leiden']\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",
    "if b_annot:\n",
    "    df_plot = df_p.groupby('leiden').mean().rename(d_replace)\n",
    "else:\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=(8,8),cmap='viridis',\n",
    "                   vmin=-2,vmax=2,method='ward')\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": [
    "ls_merge = ['20']\n",
    "#ls_merge=str_clust\n",
    "adata.obs['test'] = (adata.obs.leiden.isin(ls_merge)).replace({True:1,False:0})\n",
    "sc.pl.umap(adata, color='test',title='merged')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_merge = ['6']\n",
    "#ls_merge=str_clust\n",
    "adata.obs['test'] = (adata.obs.leiden.isin(ls_merge)).replace({True:1,False:0})\n",
    "sc.pl.umap(adata, color='test',title='merged')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot all groups spatially - leiden \n",
    "from mplex_image import analyze\n",
    "df_data = df_man\n",
    "colors = mpl.cm.tab20b.colors + mpl.cm.tab20c.colors\n",
    "\n",
    "df_pos = analyze.celltype_to_bool(df_p,'leiden')\n",
    "df_pos.columns = [str(item) for item in df_pos.columns]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_clust = '17'\n",
    "s_slide = 'RS-mTMA-5_sceneB06'#'RS-mTMA-5_sceneC03'#'RS-mTMA-5_sceneH03'#'RS-mTMA-5_sceneI10'#'cene021'\n",
    "fig,ax = plt.subplots(figsize=(5,4.5),dpi=200)\n",
    "#plot negative cells\n",
    "df_scene = df_data[df_data.index.str.contains(s_slide)]\n",
    "ax.scatter(data=df_scene,x='DAPI_X',y='DAPI_Y',color='silver',s=0.1,label=f'')\n",
    "for idxs, s_color in enumerate(df_pos.columns):\n",
    "    s_color = s_clust\n",
    "    s_color = str(s_color)\n",
    "    if len(df_p[(df_p.index.str.contains(s_slide)) & (df_pos.loc[:,s_color])])>=1:\n",
    "        #plot positive cells\n",
    "        ls_index = df_p[(df_p.index.str.contains(s_slide)) & (df_pos.loc[:,s_color])].index\n",
    "        ax.scatter(data=df_data.loc[ls_index],x='DAPI_X',y='DAPI_Y',label=f'{s_color}',s=0.1,color=colors[idxs])\n",
    "    break\n",
    "ax.set_title(f\"{s_slide}\", fontsize=16) # \\n {d_a[s_slide]}\n",
    "ax.axis('equal')\n",
    "ax.set_ylim(ax.get_ylim()[::-1])\n",
    "plt.legend(markerscale=10,framealpha=.5) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(2.5,5),dpi=200)\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_{n_neighbors}neighbors_{len(marker_genes)}markers_leiden{resolution}.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.leiden_celltype = df_p.leiden_celltype.fillna('NOS')"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## accuracy"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#merge\n",
    "df_p = df_p.merge(df_man.loc[:,['celltype']],left_index=True,right_index=True,how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.loc[df_p.celltype=='tumor','celltype'] = 'epithelial'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "#df_p['cell_type'] = df_p.celltype.fillna('NOS')\n",
    "labels = sorted(df_p.leiden_celltype.unique())\n",
    "confusion_matrix = sklearn.metrics.confusion_matrix(df_p.celltype,df_p.leiden_celltype,labels=labels)\n",
    "#accuracy\n",
    "i_acc = sklearn.metrics.accuracy_score(df_p.celltype,df_p.leiden_celltype)\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_title(f'Leiden {resolution} \\n {i_acc}')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{s_date}/Confusion_matrix_leiden{resolution}.png',dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## manual\n",
    "\n",
    "umap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "adata.obs['manual_celltype'] = df_p.celltype\n",
    "fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "figname=f'{s_sample}_celltype_manual.png'\n",
    "sc.pl.umap(adata, color='manual_celltype',ax=ax,save=figname,wspace=0.25,title=figname.split('.png')[0].replace('_',' ')) #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_date"
   ]
  },
  {
   "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",
    "#gated cells\n",
    "df_p.loc[((df_p.leiden_celltype == 'stromal') & (df_p.CD31 > 0)),'gated_celltype'] = 'endothelial'\n",
    "df_p.loc[((df_p.leiden_celltype == 'stromal') & (df_p.CD45 > 0) & (df_p.CD31 <= 0)),'gated_celltype'] = 'immune'"
   ]
  },
  {
   "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": [
    "s_sample\n",
    "resolution\n",
    "n_neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_p.merge(df_man.loc[:,['DAPI_X','DAPI_Y']],left_index=True,right_index=True).to_csv(f'results_{s_sample}_leiden{resolution}_gated_celltype.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sub Cluster  <a name=\"kmeans\"></a>\n",
    "\n",
    "re-cluster immune and fibroblast\n",
    "\n",
    "doesn't work great: gating better\n",
    " \n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_sample = '20211110_RS-mTMA-5'\n",
    "resolution = 0.5\n",
    "df = pd.read_csv(f'results_{s_sample}_leiden{resolution}_gated_celltype.csv',index_col=0)\n",
    "# cluster 5 8 and 0 are immune\n",
    "s_markers = 'stromal'\n",
    "print((df.leiden_celltype=='stromal').sum()/len(df))\n",
    "ls_marker = ['FoxP3', 'Ki67', 'CD103', 'CD11b', 'CD11c',\n",
    "       'CD31', 'CD45', 'CD8', 'ColIV', 'ColVI',  'Gal3',\n",
    "       'S100A6', 'Vim', 'aSMA', 'area', 'eccentricity', 'DAPI8', 'R1c2']\n",
    "adata = sc.AnnData(df.loc[df.leiden_celltype=='stromal',ls_marker])\n",
    "adata.raw = adata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adata.var.index\n",
    "#leiden\n",
    "for resolution in [0.6]: #0.4, 0.5, 0.6\n",
    "    results_file = f'{s_sample}_{s_markers}_markers_leiden{resolution}.h5ad'\n",
    "    print(results_file)\n",
    "    if not os.path.exists(results_file):\n",
    "        print('umap')\n",
    "        #reduce dimensionality\n",
    "        sc.tl.pca(adata, svd_solver='auto')\n",
    "        sc.pp.neighbors(adata, n_neighbors=15)\n",
    "        sc.tl.umap(adata)\n",
    "        sc.tl.leiden(adata,resolution=resolution)\n",
    "        adata.write(results_file)\n",
    "    else:\n",
    "        adata = sc.read_h5ad(results_file)\n",
    "    fig,ax = plt.subplots(figsize=(3,3),dpi=200)\n",
    "    figname=f'leiden__{s_markers}_{resolution}.png'\n",
    "    sc.pl.umap(adata, color='leiden',ax=ax,title=figname.split('.png')[0].replace('_',' ')) #save=figname,\n",
    "    #color by markers\n",
    "\n",
    "    figname = f\"Umap_subcluster_markers_{s_markers}_markers.png\"\n",
    "    title=figname.split('.png')[0].replace('_',' ')\n",
    "    sc.pl.umap(adata, color=ls_marker,vmin='p1.5',vmax='p99.5',ncols=6,save=figname)\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#clustermap\n",
    "df_p = pd.DataFrame(data=adata.X,index=adata.obs.index,columns=adata.var.index)\n",
    "df_p['leiden'] = adata.obs['leiden']\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(imm_clust,[f'{item}\\nim' for item in imm_clust])))\n",
    "#d_replace.update(dict(zip(str_clust,[f'{item}\\ns' for item in str_clust])))\n",
    "#d_replace.update(dict(zip(endo_clust,[f'{item}\\nen' for item in endo_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=(8,8),cmap='viridis',\n",
    "                   vmin=-2,vmax=2)\n",
    "#g.fig.suptitle(f'leiden {resolution}',x=.9) \n",
    "g.savefig(f'{codedir}/20210316_figures/{s_sample}_clustermap_leiden_{resolution}_{n_markers}.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": [
    "#cluster 2 and 3 are FB\n",
    "print(df.leiden.isin([2,3]).sum()/len(df))"
   ]
  },
  {
   "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 = '20211110_RS-mTMA-5'\n",
    "resolution = 0.5\n",
    "s_date = s_sample.split('_')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load annotation\n",
    "import re\n",
    "df_a = pd.read_csv('RS-MYC-PTENfl_TMA_subgroups.csv')\n",
    "df_a.loc[~df_a.GroupA.isna(),'Groups'] = 'stroma_rich'\n",
    "df_a.loc[~df_a.GroupB.isna(),'Groups'] = 'stroma_poor'\n",
    "df_a.loc[df_a.Groups.isna(),'Groups'] = 'not_tumor'\n",
    "df_a['Xiaoyan'] = df_a.Xiaoyan.fillna(df_a.TMA)\n",
    "df_a.index = df_a.TMA #drop(['oldE1','oldF6'],inplace=True)\n",
    "d_a = dict(zip(df_a.index,df_a.Groups))\n",
    "df_a['slide_scene'] = [f'RS-mTMA-5_scene{re.sub(r\"[0-9]+\", \"\", item)}{int(re.sub(r\"[A-Z]+\", \"\", item)):02d}' for item in df_a.index]\n",
    "df_a['Group'] = df_a.Xiaoyan.map(d_a).fillna('not_tumor')      \n",
    "\n"
   ]
  },
  {
   "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]\n",
    "#annotate\n",
    "d_a = dict(zip(df_a.slide_scene,df_a.Group))\n",
    "df_p['Group'] = df_p.slide_scene.map(d_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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_scene)):\n",
    "    s_slide = 'RS-mTMA-5_sceneB07'\n",
    "    s_slide = 'RS-mTMA-5_sceneA09'\n",
    "    fig,ax = plt.subplots(figsize=(5,4.5),dpi=200)\n",
    "    #plot negative cells\n",
    "    df_scene = df_p[df_p.index.str.contains(s_slide)]\n",
    "    ax.scatter(data=df_scene,x='DAPI_X',y='DAPI_Y',color='silver',s=0.1,label=f'')\n",
    "    #for idxs, s_color_int in enumerate(range(len(df_pos.columns))):\n",
    "    for idxs, s_color in enumerate(df_pos.columns):\n",
    "        s_color = str(s_color)\n",
    "        if len(df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ax.scatter(data=df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])],x='DAPI_X',y='DAPI_Y',\n",
    "                                                                        label=f'{s_color}',s=0.1,color=mpl.cm.tab20.colors[idxs])\n",
    "        #break\n",
    "    ax.set_title(f\"{s_slide} \\n {d_a[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) \n",
    "    fig.savefig(f'{codedir}/{s_date}/{s_slide}_leiden{resolution}_scatterplot.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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_scene)):\n",
    "    s_slide = 'RS-mTMA-5_sceneB07'\n",
    "    s_slide = 'RS-mTMA-5_sceneA09'\n",
    "    fig,ax = plt.subplots(figsize=(5,5),dpi=300)\n",
    "    #plot negative cells\n",
    "    df_scene = df_p[df_p.index.str.contains(s_slide)]\n",
    "    ax.scatter(data=df_scene,x='DAPI_X',y='DAPI_Y',color='silver',s=0.1,label=f'')\n",
    "    #for idxs, s_color_int in enumerate(range(len(df_pos.columns))):\n",
    "    for idxs, s_color in enumerate(df_pos.columns):\n",
    "        #print(s_color)\n",
    "        #s_color = str(s_color_int)\n",
    "        if len(df_p[(df_p.slide_scene==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ax.scatter(data=df_p[(df_p.slide_scene==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} \\n {d_a[s_slide]}\", fontsize=16)\n",
    "    ax.axis('equal')\n",
    "    if s_slide == 'RS-mTMA-5_sceneB07':\n",
    "        ax.set_ylim(400,5350)\n",
    "        ax.set_xlim(400,5350)\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    \n",
    "    #ax.set_xticklabels('')\n",
    "    #ax.set_yticklabels('')\n",
    "    #break\n",
    "    plt.legend(markerscale=10,framealpha=.1) \n",
    "    fig.savefig(f'{codedir}/{s_date}/{s_slide}_gated_cells_scatterplot.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## subtype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_type = 'leiden' #'gated_celltype' #\n",
    "df = df_p.loc[:,['CD31',s_type,'slide_scene']].groupby([s_type,'slide_scene']).CD31.count()/df_p.loc[:,['CD31','slide_scene']].groupby(['slide_scene']).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,size=300,save=figname)\n",
    "    #survival\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",
    "   \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": {},
   "outputs": [],
   "source": [
    "#sns.clustermap(df_prop)\n",
    "ls_drop = df_prop[df_prop.index.isin(df_a[df_a.Group=='not_tumor'].slide_scene)].index\n",
    "df_prop_plot = df_prop.drop(ls_drop)\n",
    "\n",
    "\n",
    "d_color = dict(zip(['stroma_rich','stroma_poor'],sns.color_palette()[0:2]))\n",
    "row_colors = df_prop_plot.index.map(d_a).map(d_color)\n",
    "\n",
    "\n",
    "g = sns.clustermap(df_prop_plot,figsize=(6,6),cmap='viridis',row_colors=row_colors,\n",
    "        )\n",
    "#g.fig.suptitle(f'leiden {resolution}',x=.9) \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}_{s_type}_clustermap_Group_leiden{resolution}.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()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prop.loc[:,df_p.loc[df_p.celltype=='epithelial','leiden'].unique()]\n",
    "df_p.loc[df_p.celltype=='epithelial','leiden'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_prop_plot[df_prop_plot.sum(axis=1)>0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_cell in ['stromal','epithelial','immune']: #'immune'\n",
    "    ls_index = df_prop[~df_prop.index.isin(df_a[df_a.Group=='not_tumor'].slide_scene)].index\n",
    "    ls_col = df_p.loc[df_p.gated_celltype==s_cell,'leiden'].unique()\n",
    "    df_prop_plot = df_prop.loc[ls_index,ls_col]\n",
    "    df_prop_plot = df_prop_plot[df_prop_plot.sum(axis=1)>0]\n",
    "    d_color = dict(zip(['stroma_rich','stroma_poor'],sns.color_palette()[0:2]))\n",
    "    row_colors = df_prop_plot.index.map(d_a).map(d_color)\n",
    "\n",
    "    g = sns.clustermap(df_prop_plot,figsize=(6,6),cmap='viridis',row_colors=row_colors,z_score=0)\n",
    "            #vmin=-2,vmax=2)\n",
    "    g.fig.suptitle(f'{s_cell}',x=.7) \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}_{s_type}_clustermap_Group_leiden{resolution}_{s_cell}.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": [
    "#leiden\n",
    "for resolution in [0.1]: #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": {},
   "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=(3,4),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": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}