Normalize_Bx2-4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load libraries\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import os\n",
    "import copy\n",
    "import seaborn as sns\n",
    "import importlib\n",
    "from scipy.signal import argrelmax, find_peaks, peak_widths\n",
    "import scanpy as sc\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.preprocessing import scale, minmax_scale\n",
    "from sklearn.metrics import silhouette_score\n",
    "import matplotlib as mpl\n",
    "mpl.rc('figure', max_open_warning = 0)\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "mpl.rcParams['mathtext.it'] = 'Arial:italic'\n",
    "mpl.rcParams['mathtext.rm'] = 'Arial'\n",
    "codedir = os.getcwd()\n",
    "#load cmif libraries\n",
    "#os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cmIF')\n",
    "from mplex_image import visualize as viz, process, preprocess, normalize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(codedir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(222)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of contents <a name=\"contents\"></a>\n",
    "1. [Load Data](#load)\n",
    "2. [Normalize](#norm)\n",
    "3. [Visualize Normalization](#normviz)\n",
    "4. [leiden for cell typing](#clusterlei)\n",
    "5. [Leiden cluster](#clust1)\n",
    "\n",
    "\n",
    "note:\n",
    "\n",
    "    Could you make composite fraction bar graph only  in following regions?\n",
    "\n",
    "    Bx2: SMTBx2-5-Scene-001_ROI1-2000-9000-2500-2500\n",
    "    Bx3: SMTBx3-Scene-004_ROI2-20900-15494-2500-2500\n",
    "    Bx4: HTA-33-Scene-002_ROI1-3271-607-2500-2500\n",
    "\n",
    "    If we can have it in Bx1\n",
    "    Bx: SMTBx1-Scene-003_ROI1-2440-220-2500-2500\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data\n",
    "os.chdir(f'{codedir}/paper_data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_date = '20210402'\n",
    "if not os.path.exists(s_date):\n",
    "    os.mkdir(s_date)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Data <a name=\"load\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(f'{codedir}/paper_data')\n",
    "df_file = pd.DataFrame(index=os.listdir())\n",
    "df_file = df_file[df_file.index.str.contains('FilteredMeanIntensity_DAPI')]\n",
    "df_file['tissue'] = [item.split('_')[1] for item in df_file.index]\n",
    "df_file['dapi'] = ['DAPI' + item.split('y_DAPI')[1].split('.')[0] for item in df_file.index]\n",
    "ls_sample = df_file.tissue.tolist()\n",
    "d_dapi = dict(zip(df_file.tissue.tolist(),df_file.dapi.tolist()))\n",
    "d_dapi.update({'JE-TMA-60': 'DAPI10_DAPI2'})\n",
    "df_mi = pd.DataFrame()\n",
    "df_xy = pd.DataFrame()\n",
    "df_edge = pd.DataFrame()\n",
    "\n",
    "for s_sample in sorted(set(ls_sample)):\n",
    "    #if not s_sample.find('HTA')>-1:\n",
    "        print(f'loading {s_sample}')\n",
    "        df_mi = df_mi.append(pd.read_csv(f'{codedir}/paper_data/features_{s_sample}_FilteredMeanIntensity_{d_dapi[s_sample]}.csv', index_col=0))\n",
    "        df_xy = df_xy.append(pd.read_csv(f'{codedir}/paper_data/features_{s_sample}_CentroidXY.csv',index_col=0))\n",
    "        if os.path.exists(f'{codedir}/paper_data/features_{s_sample}_EdgeCells153pixels_CentroidXY.csv'):\n",
    "            df_edge = df_edge.append(pd.read_csv(f'{codedir}/paper_data/features_{s_sample}_EdgeCells153pixels_CentroidXY.csv',index_col=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sorted(df_mi.columns[df_mi[~df_mi.index.str.contains('JE-TMA-60')].isna().sum() != 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_marker = ['AR_nuclei', 'CD20_perinuc5', 'CD31_perinuc5', 'CD3_perinuc5',  'CD44_perinuc5', 'CD45_perinuc5',#'CD44_nucadj2',\n",
    " 'CD4_perinuc5', 'CD68_perinuc5','CD8_perinuc5', 'CK14_cytoplasm', 'CK17_cytoplasm', 'CK19_cytoplasm', 'CK5_cytoplasm',\n",
    " 'CK7_cytoplasm', 'CK8_cytoplasm', 'ColI_perinuc5', 'ColIV_perinuc5','CoxIV_perinuc5','EGFR_cytoplasm', 'ER_nuclei',\n",
    " 'Ecad_cytoplasm', 'FoxP3_nuclei', 'GRNZB_nuclei', 'H3K27_nuclei','H3K4_nuclei', 'HER2_cellmem25','Ki67_nuclei',\n",
    " 'LamAC_nuclei', 'PCNA_nuclei', 'PD1_perinuc5', 'PDPN_perinuc5','DAPI2_nuclei',  # 'ER_nuclei25','HER2_cytoplasm','PgR_nuclei','Vim_nucadj2'\n",
    " 'Vim_perinuc5', 'aSMA_perinuc5', 'pHH3_nuclei', 'pRB_nuclei', 'pS6RP_perinuc5','slide_scene',\n",
    "           ] # CD8R bad, 'gH2AX_nuclei' in R11 Bx3 not included\n",
    "\n",
    "df_mi = df_mi.loc[:,ls_marker]\n",
    " \n",
    "# old \n",
    "#df_mi = df_mi.loc[:,['HER2_cellmem25', 'DAPI2_nuclei',# 'CD44_nucadj2', 'Vim_nucadj2','ER_nuclei25','HER2_cytoplasm',\n",
    "#       'CD20_perinuc5', 'CD3_perinuc5', 'CD31_perinuc5', 'CD4_perinuc5','CD44_perinuc5', 'CD45_perinuc5', 'CD68_perinuc5', 'CD8_perinuc5',\n",
    "#       'PD1_perinuc5', 'PDPN_perinuc5', 'Vim_perinuc5', 'aSMA_perinuc5','CK14_cytoplasm', 'CK17_cytoplasm', 'CK19_cytoplasm', 'CK5_cytoplasm',\n",
    "#       'CK7_cytoplasm', 'Ecad_cytoplasm', 'ER_nuclei', 'Ki67_nuclei', 'LamAC_nuclei','PCNA_nuclei', 'pHH3_nuclei', 'slide_scene']]\n",
    "\n",
    "\n",
    "df_mi['batch'] = [item.split('_')[0] for item in df_mi.index]\n",
    "#df_mi['scene'] = [item.split('_')[1] for item in df_mi.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deal with JE-TMA-60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# markers in JE-TMA-60\n",
    "#'JE-TMA-60_scene06', 'JE-TMA-60_scene08', 'JE-TMA-60_scene09', 'JE-TMA-60_scene10', 'JE-TMA-60_scene11', 'JE-TMA-60_scene13'\n",
    "# R5 is CK17.PDPN.CD45.FoxP3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_R5 = pd.read_csv(f'{codedir}/paper_data/features_JE-TMA-60_FilteredMeanIntensity_DAPI5_DAPI2.csv',index_col=0)\n",
    "df_R4 = pd.read_csv(f'{codedir}/paper_data/features_JE-TMA-60_FilteredMeanIntensity_DAPI4_DAPI2.csv',index_col=0)\n",
    "df_R10 = df_mi[df_mi.batch=='JE-TMA-60']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_scene = set(df_R10.slide_scene)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_na = set([item.split('_cell')[0] for item in df_R5.index]) - set([item.split('_cell')[0] for item in df_R10.index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#slect markers, scenes for normalization (based on JE-TMA-60 tissue loss)\n",
    "ls_pos = ['HER2_cellmem25','CK19_cytoplasm','CK7_cytoplasm','CK8_cytoplasm','Ecad_cytoplasm','ER_nuclei','Ki67_nuclei','LamAC_nuclei',\n",
    "          'PCNA_nuclei','pHH3_nuclei','Vim_perinuc5','DAPI2_nuclei','H3K27_nuclei','H3K4_nuclei', 'pRB_nuclei','pS6RP_perinuc5',\n",
    "         'CoxIV_perinuc5','EGFR_cytoplasm']\n",
    "ls_R5 = ['CK17_cytoplasm','PDPN_perinuc5','CD45_perinuc5','FoxP3_nuclei'] #\n",
    "ls_R4 = ['pHH3_nuclei','CK14_cytoplasm','Ki67_nuclei','CK19_cytoplasm','CK5_cytoplasm','HER2_cellmem25',\n",
    "        'Ecad_cytoplasm', 'ER_nuclei','CD44_perinuc5', 'PCNA_nuclei','aSMA_perinuc5','CD3_perinuc5','EGFR_cytoplasm']\n",
    "ls_bad = ['CD20_perinuc5', 'CD31_perinuc5', 'CD4_perinuc5', 'CD68_perinuc5', 'CD8_perinuc5','PD1_perinuc5',\n",
    "         'ColI_perinuc5', 'ColIV_perinuc5'] #'CK7_cytoplasm', #'LamAC_nuclei',\n",
    "#ls_good = ['CK7_cytoplasm','Vim_perinuc5','LamAC_nuclei']\n",
    "\n",
    "#R4\n",
    "df =  df_mi[df_mi.batch!='JE-TMA-60']\n",
    "df = df.append(df_R4.loc[:,ls_R4])\n",
    "#R5\n",
    "ls_index = df_R5.loc[df_R5.index.isin(df_R4.index)].index\n",
    "df.loc[ls_index,ls_R5] = df_R5.loc[ls_index,ls_R5]\n",
    "\n",
    "#fill R6-8\n",
    "ls_index = df_mi.loc[(df_mi.slide_scene.isin(ls_scene)) & (df_mi.index.isin(df_R4.index))].index\n",
    "df.loc[ls_index,ls_pos] = df_R10.loc[ls_index,ls_pos]\n",
    "\n",
    "#\n",
    "df['batch'] = [item.split('_')[0] for item in df.index]\n",
    "#df['scene'] = [item.split('_')[1] for item in df.index]\n",
    "df['slide_scene'] = [item.split('_cell')[0] for item in df.index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## filter edge cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#filter out unwanted cells\n",
    "d_filter = {#41 (not used)\n",
    "            'JE-TMA-41_scene01':(df_xy.DAPI_Y > 5000),'JE-TMA-41_scene03':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-41_scene04':(df_xy.DAPI_Y < 1500),'JE-TMA-41_scene05':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-41_scene06':(df_xy.DAPI_Y < 1500),'JE-TMA-41_scene08':(df_xy.DAPI_Y < 1500),\n",
    "            'JE-TMA-41_scene09':(df_xy.DAPI_Y > 5000),'JE-TMA-41_scene11':(df_xy.DAPI_Y < 1500),\n",
    "            #43\n",
    "            'JE-TMA-43_scene09':(df_xy.DAPI_Y < 1200),'JE-TMA-43_scene14':(df_xy.DAPI_Y < 1200),\n",
    "            #60\n",
    "            'JE-TMA-60_scene02':(df_xy.DAPI_X < 1500),'JE-TMA-60_scene05':(df_xy.DAPI_X < 1500),\n",
    "            'JE-TMA-60_scene11':(df_xy.DAPI_Y < 1500),'JE-TMA-60_scene14':(df_xy.DAPI_X < 1500),\n",
    "            'JE-TMA-60_scene06':(df_xy.DAPI_Y < 1500),'JE-TMA-60_scene08':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-60_scene10':(df_xy.DAPI_Y < 1500),\n",
    "            #63\n",
    "            'JE-TMA-62_scene01':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-62_scene02':(df_xy.DAPI_X > 5000),'JE-TMA-62_scene03':(df_xy.DAPI_X < 1000),\n",
    "            'JE-TMA-62_scene04':(df_xy.DAPI_Y < 1500),'JE-TMA-62_scene06':(df_xy.DAPI_X < 1000),\n",
    "            'JE-TMA-62_scene08':(df_xy.DAPI_Y > 5000),'JE-TMA-62_scene10':(df_xy.DAPI_Y < 1500),\n",
    "            #'SMTBx1-16_scene001':(df_xy.DAPI_Y > 1), #keep scene 1 for manual thresholding\n",
    "            'SMTBx2-3_scene002':(df_xy.DAPI_Y > 5000),'SMTBx3_scene004':(df_xy.DAPI_X <11000),\n",
    "            'SMTBx3_scene005':(df_xy.DAPI_X > 0),'SMTBx4-3_scene001':(df_xy.DAPI_Y < 2400),\n",
    "            'SMTBx2-5_scene002':(df_xy.DAPI_Y > 5000),'HTA-33_scene003':(df_xy.DAPI_Y > 9000)}\n",
    "d_filter2 = {'JE-TMA-60_scene02':(df_xy.DAPI_Y > 4500)}\n",
    "ls_filter_all = []\n",
    "for s_scene, filtercon in d_filter.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "for s_scene, filtercon in d_filter2.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "#filter edge\n",
    "ls_filter_all = ls_filter_all + df_edge.index.tolist()\n",
    "df_filter_mi = df[(~df.index.isin(ls_filter_all))]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_cluster = df_filter_mi.loc[:,['HER2_cellmem25','slide_scene']]\n",
    "df_cluster['cluster'] = 1\n",
    "df_cluster.drop('HER2_cellmem25',axis=1,inplace=True)\n",
    "import importlib\n",
    "importlib.reload(viz)\n",
    "%matplotlib inline\n",
    "viz.plot_clusters(df_cluster,df_xy,s_num='few')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#match controls to biopsies\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "d_replace = {'BC44290-146': 'JE-TMA-41',\n",
    " 'SMTBx2-3': 'JE-TMA-41',\n",
    " 'SMTBx2-5':'JE-TMA-43',\n",
    " 'SMTBx3':'JE-TMA-60',\n",
    " 'SMTBx4-3':'JE-TMA-62'}\n",
    "df_filter_mi.loc[:,'batch'] = df_filter_mi.batch.replace(d_replace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#standardize the scenes\n",
    "d_replace = { 'JE-TMA-41_scene13':'JE-TMA-41_scene14',\n",
    "             'JE-TMA-41_scene12':'JE-TMA-41_scene13',\n",
    "             'JE-TMA-62_scene13':'JE-TMA-62_scene14',\n",
    "             'JE-TMA-62_scene12':'JE-TMA-62_scene13'}\n",
    "df_filter_mi.loc[:,'scene'] = df_filter_mi.slide_scene.replace(d_replace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_filter_mi.merge(df_xy.loc[:,['DAPI_X', 'DAPI_Y', 'nuclei_area', 'nuclei_eccentricity']],left_index=True,right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_out = df_filter_mi.merge(df_xy.loc[:,['DAPI_X', 'DAPI_Y', 'nuclei_area', 'nuclei_eccentricity']],left_index=True,right_index=True)\n",
    "len(df_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2-23 contains NAs\n",
    "#2-22 the NAs were filled with random gaussian data\n",
    "# 0302 include scene 1 Bx1\n",
    "# 0318 just Bx2 - 4, (Bx2-5)\n",
    "# 20210324 has HTA9-1-33\n",
    "if not os.path.exists('20210324_SMTBx1-4_JE-TMA-43_60_62_FilteredMeanIntensity.csv'):\n",
    "    print('saving csv')\n",
    "    #df_out.to_csv('20210223_SMTBx1-4_JE-TMA-41_60_62_BC44290-146.csv')\n",
    "    df_out.to_csv('20210324_SMTBx1-4_JE-TMA-43_60_62_FilteredMeanIntensity.csv') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2-23 contains NAs\n",
    "#2-22 the NAs were filled with random gaussian data\n",
    "# 0302 include scene 1 Bx1\n",
    "# 0318 just Bx2 - 4, (Bx2-5)\n",
    "if not os.path.exists('20210320_SMTBx2-4_JE-TMA-43_60_62_FilteredMeanIntensity.csv'):\n",
    "    print('saving csv')\n",
    "    #df_out.to_csv('20210223_SMTBx1-4_JE-TMA-41_60_62_BC44290-146.csv')\n",
    "    df_out.to_csv('20210320_SMTBx2-4_JE-TMA-43_60_62_FilteredMeanIntensity.csv') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization <a name=\"norm\"></a>\n",
    "\n",
    "use ComBat.\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi = pd.read_csv('20210320_SMTBx2-4_JE-TMA-43_60_62_FilteredMeanIntensity.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mi.scene.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_pos = ['HER2_cellmem25','CK19_cytoplasm','CK7_cytoplasm','CK8_cytoplasm','Ecad_cytoplasm','ER_nuclei','Ki67_nuclei','LamAC_nuclei',\n",
    "          'PCNA_nuclei','pHH3_nuclei','Vim_perinuc5','DAPI2_nuclei','H3K27_nuclei','H3K4_nuclei', 'pRB_nuclei','pS6RP_perinuc5',\n",
    "         'CoxIV_perinuc5','EGFR_cytoplasm']\n",
    "ls_R5 = ['CK17_cytoplasm','PDPN_perinuc5','CD45_perinuc5','FoxP3_nuclei'] #\n",
    "ls_R4 = ['pHH3_nuclei','CK14_cytoplasm','Ki67_nuclei','CK19_cytoplasm','CK5_cytoplasm','HER2_cellmem25',\n",
    "        'Ecad_cytoplasm', 'ER_nuclei','CD44_perinuc5', 'PCNA_nuclei','aSMA_perinuc5','CD3_perinuc5','EGFR_cytoplasm']\n",
    "ls_bad = ['CD20_perinuc5', 'CD31_perinuc5', 'CD4_perinuc5', 'CD68_perinuc5', 'CD8_perinuc5','PD1_perinuc5',\n",
    "         'ColI_perinuc5', 'ColIV_perinuc5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select normalization scenes\n",
    "ls_R10_scene = ['scene06', 'scene08', 'scene09', 'scene10', 'scene11', 'scene13']\n",
    "ls_R10 = ['HER2_cellmem25', 'CK19_cytoplasm', 'CK7_cytoplasm', 'Ecad_cytoplasm', 'ER_nuclei', 'Ki67_nuclei', 'LamAC_nuclei',\n",
    "          'PCNA_nuclei','pHH3_nuclei', 'Vim_perinuc5','CD44_perinuc5','DAPI2_nuclei', #adding following:\n",
    "          'CK8_cytoplasm','CoxIV_perinuc5', 'EGFR_cytoplasm', 'H3K27_nuclei', 'H3K4_nuclei', 'pRB_nuclei', 'pS6RP_perinuc5']\n",
    "#note: CK17 may have quenching artifact; PDPN not good in Bx1, so just CD45 important\n",
    "#'CK17_cytoplasm','PDPN_perinuc5', #'FoxP3_nuclei' not in full set\n",
    "ls_R5 = ['PDPN_perinuc5','CD45_perinuc5','FoxP3_nuclei', 'aSMA_perinuc5','CD3_perinuc5']  # aSMA because N breast, scene 01 better than 07 for immune\n",
    "ls_R5_scene = ['scene01','scene03','scene04']\n",
    "#old ls_R4 = ['pHH3_nuclei','CK14_cytoplasm','Ki67_nuclei','CK19_cytoplasm','CK5_cytoplasm','HER2_cellmem25',\n",
    "#        'Ecad_cytoplasm', 'ER_nuclei','CD44_perinuc5', 'PCNA_nuclei','aSMA_perinuc5','CD3_perinuc5','DAPI2_nuclei']\n",
    "#can scene 7 be good control for CD3 and CK14 and CK5?, yes. R1 doen't add much\n",
    "ls_R4 = [ 'CK14_cytoplasm', 'CK5_cytoplasm','CK17_cytoplasm'] #'CD3_perinuc5',\n",
    "ls_R4_scene = ['scene02','scene07']\n",
    "ls_bad = ['CD20_perinuc5', 'CD31_perinuc5', 'CD4_perinuc5', 'CD68_perinuc5', 'CD8_perinuc5','PD1_perinuc5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set(df_mi.batch)\n",
    "#df_mi = df_mi.loc[df_mi.batch!='JE-TMA-60']\n",
    "df_mi['slide_scene'] = df_mi.scene\n",
    "df_mi['scene'] = [item.split('_')[1] for item in df_mi.slide_scene]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dropped 60\n",
    "df_norm_all=pd.DataFrame(index=df_mi.dropna().index)\n",
    "\n",
    "#not dropped 60\n",
    "df_norm_all=pd.DataFrame(index=df_mi.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1 fit on scenes that are good through round 10 and markers that are positive on those scenes \"pos\"\n",
    "for s_type in ['R4','R5','R10']:\n",
    "    if s_type == 'R10':\n",
    "        ls_pos = ls_R10\n",
    "        ls_scene = ls_R10_scene\n",
    "\n",
    "    #2 fit on scenes that are good until R4, and R1-4 markers\n",
    "    if s_type == 'R4':\n",
    "        ls_pos = ls_R4\n",
    "        ls_scene = ls_R4_scene # + ls_R5_scene + ls_R10_scene \n",
    "\n",
    "    #3 fit on scene that are good until R5, and R5 markers\n",
    "    if s_type == 'R5':\n",
    "        ls_pos = ls_R5\n",
    "        ls_scene = ls_R5_scene\n",
    "\n",
    "    #fit\n",
    "    b_control = ((df_mi.index.str.contains('JE-TMA')) & (df_mi.scene.isin(ls_scene)) & (df_mi.loc[:,ls_pos].isna().sum(axis=1)==0))\n",
    "    data = df_mi.loc[b_control,ls_pos].T\n",
    "    batch = df_mi.loc[b_control,'batch']\n",
    "    gamma_star, delta_star, stand_mean, var_pooled = normalize.combat_fit(data, batch)\n",
    "    #transform\n",
    "    #data = df_mi.loc[df_mi.batch!='SMTBx1-16',df_mi.dtypes=='float64'].drop(['DAPI_X','DAPI_Y'],axis=1).T\n",
    "    data = df_mi.loc[df_mi.batch!='SMTBx1-16',ls_pos].T\n",
    "    batch = df_mi.loc[df_mi.batch!='SMTBx1-16','batch']\n",
    "    bayesdata = normalize.combat_transform(data,batch,gamma_star,delta_star,stand_mean, var_pooled)\n",
    "    df_norm = bayesdata.T\n",
    "    df_norm_all = df_norm_all.merge(df_norm,left_index=True,right_index=True,how='left')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_norm_all.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run after #1, 2 and 3\n",
    "df_norm_all =  df_norm_all.merge(df_mi.loc[:,['batch','DAPI_X','DAPI_Y','scene','nuclei_area','nuclei_eccentricity']],left_index=True,right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#old check\n",
    "df_norm = df_norm.merge(df_mi.loc[:,['batch','DAPI_X','DAPI_Y','scene','nuclei_area','nuclei_eccentricity']],left_index=True,right_index=True)\n",
    "#df_mi.loc[b_control,:].drop(['DAPI_X','DAPI_Y'],axis=1).groupby('batch').mean()\n",
    "#df_mi[df_mi.index.str.contains('JE-TMA')].drop(['DAPI_X','DAPI_Y'],axis=1).groupby('batch').std()\n",
    "#check\n",
    "df_norm.loc[b_control,:].drop(['DAPI_X','DAPI_Y'],axis=1).groupby('batch').mean()\n",
    "#df_norm[df_norm.index.str.contains('JE-TMA')].drop(['DAPI_X','DAPI_Y'],axis=1).groupby('batch').std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_norm_all.to_csv('20210320_SMTBx2-4_JE-TMA-43_60_62_normalized.csv')\n",
    "#df_norm_all.to_csv('20210325_SMTBx2-4_JE-TMA-43_60_62_normalized.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Umap Visualize Normalization <a name=\"normviz\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#s_sample = '20210320_SMTBx2-4_JE-TMA-43_60_62'\n",
    "s_sample = '20210325_SMTBx2-4_JE-TMA-43_60_62'\n",
    "df_norm_all = pd.read_csv(f'{s_sample}_normalized.csv',index_col=0)\n",
    "df_norm_all.rename({'nuclei_area':'area','nuclei_eccentricity':'eccentricity','DAPI_X':'DAPIX',\n",
    "       'DAPI_Y':\"DAPIY\"},axis=1, inplace=True)\n",
    "df_norm_all.columns = [item.split('_')[0] for item in df_norm_all.columns]\n",
    "df_norm_all['slide'] = [item.split('_')[0] for item in df_norm_all.index]\n",
    "df_norm_all['scene'] = [item.split('_')[1] for item in df_norm_all.index]\n",
    "df_norm_all['slide_scene'] = [item.split('_cell')[0] for item in df_norm_all.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_norm_all = df_norm_all.loc[~df_norm_all.slide_scene.isin(['JE-TMA-43_scene01','JE-TMA-62_scene01'])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# visualize\n",
    "%matplotlib inline\n",
    "s_type = 'w-60_no01'\n",
    "#adata = sc.AnnData(df_norm_all.loc[:,df_norm_all.dtypes=='float64'].drop(['DAPIX','DAPIY'],axis=1)) \n",
    "ls_drop = ['DAPIX','DAPIY','DAPI2','LamAC','pHH3','FoxP3','CoxIV',\n",
    "          'H3K27','H3K4','pRB','pS6RP','aSMA','PDPN'] #aSMA, PDPN not well norm\n",
    "adata = sc.AnnData(df_norm_all.dropna().loc[:,df_norm_all.dtypes=='float64'].drop(ls_drop,axis=1))\n",
    "adata.obs['batch'] = df_norm_all.dropna().loc[:,'batch']\n",
    "adata.obs['scene'] =  df_norm_all.dropna().loc[:,'scene'].replace({'scene001':'Bx', 'scene002':'Bx','scene003':'Bx', 'scene004':'Bx', 'scene005':'Bx'})\n",
    "adata.obs['tissue'] = df_norm_all.dropna().loc[:,'slide']\n",
    "# reduce dimensionality (PCA)\n",
    "adata.raw = adata\n",
    "#reduce dimensionality\n",
    "sc.tl.pca(adata, svd_solver='auto')\n",
    "#sc.pl.pca(adata)\n",
    "sc.pl.pca_variance_ratio(adata, log=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate neighbors\n",
    "n_neighbors = 31\n",
    "n_pcs=len(adata.var.index) - 1\n",
    "results_file =  f'{s_sample}_{n_neighbors}neighbors_{n_pcs}pcs_{len(adata.var.index)}markers.h5ad'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_celline = {'scene02':'HCC1143',\n",
    " 'scene03':'HCC3153',\n",
    " 'scene04':'N.Breast',\n",
    " 'scene05':'T47D',\n",
    " 'scene06':'T47D',\n",
    " 'scene07':'Tonsil',\n",
    " 'scene08':'BT474',\n",
    " 'scene09':'BT474',\n",
    " 'scene10':'AU565',\n",
    " 'scene11':'AU565',\n",
    " 'scene12':'MDAMB436',\n",
    " 'scene13':'MDAMB436',\n",
    " 'scene14':'MDAMB436'}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# calculate neighbors\n",
    "if os.path.exists(results_file):\n",
    "    adata = sc.read_h5ad(results_file)\n",
    "    print('loading umap')\n",
    "else:\n",
    "    # calculate neighbors \n",
    "    print('calculating umap')\n",
    "    sc.pp.neighbors(adata, n_neighbors=n_neighbors, n_pcs=n_pcs)\n",
    "    sc.tl.umap(adata)\n",
    "    #save results\n",
    "    if not os.path.exists(results_file):\n",
    "        adata.write(results_file)\n",
    "\n",
    "# umap plus scenes\n",
    "fig,ax = plt.subplots(figsize=(3,2.5),dpi=600)\n",
    "figname = f'UmapScene_{s_type}_{n_pcs+1}markers.png'\n",
    "sc.pl.umap(adata, color='scene',save=figname,title=f'TMA Core',ax=ax)\n",
    "\n",
    "\n",
    "# umap plus tissue\n",
    "fig,ax = plt.subplots(figsize=(3,2.5),dpi=600)\n",
    "figname = f'UmapTissue_{s_type}_{n_pcs+1}markers.png'\n",
    "adata.obs['Tissue'] = adata.obs['tissue'].replace({'SMTBx2-5':'Bx2', 'SMTBx3':'Bx3','SMTBx4-3':'Bx4'})\n",
    "sc.pl.umap(adata, color='Tissue',save=figname,title=f'Tissue',ax=ax)\n",
    "\n",
    "\n",
    "# umap plus cell line\n",
    "adata.obs['Subtype'] = adata.obs.scene.replace(d_celline)\n",
    "fig,ax = plt.subplots(figsize=(3,2.5),dpi=600)\n",
    "figname = f'UmapSubtype_{s_type}_{n_pcs+1}markers.png'\n",
    "sc.pl.umap(adata, color='Subtype',save=figname,title=f'Subtype',ax=ax)\n",
    "\n",
    "\n",
    "#umap plot\n",
    "ls_marker = adata.var.index.tolist()\n",
    "figname = f\"Umap_{s_type}_{n_pcs+1}markers.png\"\n",
    "axes = sc.pl.umap(adata, color=ls_marker,wspace=.25,save=figname,vmin='p1.5',vmax='p98.5',ncols=3,show=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#umap plot\n",
    "ls_marker = adata.var.index.tolist()\n",
    "figname = f\"Umap_{s_type}_{n_pcs+1}markers.png\"\n",
    "fig = sc.pl.umap(adata, color=ls_marker,wspace=.25,vmin='p1.5',vmax='p98.5',ncols=3,show=False,return_fig=True)\n",
    "ax_list = fig.axes\n",
    "for ax in ax_list:\n",
    "    ax.set_title(ax.get_title(),fontsize=28)\n",
    "fig.savefig(f'figures/{figname}',dpi=600)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## cluster leiden <a name=\"clusterlei\"></a>\n",
    "\n",
    "[contents](#contents)\n",
    "\n",
    "cluster on the markers that are normalized well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "resolution = 0.45\n",
    "results_file = f'{s_sample}_{n_neighbors}neighbors_{n_pcs}pcs_{len(adata.var.index)}markers_leiden{resolution}.h5ad'\n",
    "#save\n",
    "if not os.path.exists(results_file):\n",
    "    sc.tl.leiden(adata,resolution=resolution)\n",
    "else:\n",
    "    adata = sc.read_h5ad(results_file)\n",
    "    print('loading leiden')    \n",
    "fig,ax = plt.subplots(figsize=(3,2.5),dpi=600)\n",
    "figname=f'leiden_{resolution}.png'\n",
    "sc.pl.umap(adata, color='leiden',ax=ax,save=figname)\n",
    "#seaborn 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",
    "g = sns.clustermap(df_p.groupby('leiden').mean(),z_score=1,figsize=(4,4),cmap='viridis',\n",
    "                   vmin=-1.5,vmax=1.5) \n",
    "#g.savefig(f'./figures/clustermap_leiden.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()\n",
    "#scanpy matrixplot\n",
    "fig,ax = plt.subplots(figsize=(5,5), dpi=200)\n",
    "figname=f'Matrixplot_leiden_{resolution}.png'\n",
    "sc.pl.matrixplot(adata, var_names=marker_genes, groupby=f'leiden',title='',categories_order=categories_order,\n",
    "                 ax=ax,save=figname,standard_scale='var',colorbar_title='Relative\\nintensity',\n",
    "                #var_group_positions=[(3,23),(24,31),(32,42),(43,51)],\n",
    "                 #var_group_labels=['tumor','T-cell','muscle\\n +AF','immune\\n+stroma'],\n",
    "                #var_group_rotation=0\n",
    "                )\n",
    "\n",
    "#save\n",
    "if not os.path.exists(results_file):\n",
    "    adata.write(results_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leiden barplots <a name=\"clust1\"></a>\n",
    "\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_order = [\n",
    "       'Bx2','Bx3','Bx4',#'JE-TMA-43_scene01','JE-TMA-62_scene01',\n",
    "       'JE-TMA-43_scene02', 'JE-TMA-62_scene02',\n",
    "       'JE-TMA-43_scene03', 'JE-TMA-62_scene03', 'JE-TMA-43_scene04',\n",
    "       'JE-TMA-62_scene04', 'JE-TMA-43_scene05', 'JE-TMA-62_scene05',\n",
    "       'JE-TMA-43_scene06','JE-TMA-60_scene06', 'JE-TMA-62_scene06', 'JE-TMA-43_scene07',\n",
    "       'JE-TMA-62_scene07', 'JE-TMA-43_scene08','JE-TMA-60_scene08', 'JE-TMA-62_scene08',\n",
    "       'JE-TMA-43_scene09','JE-TMA-60_scene09', 'JE-TMA-62_scene09','JE-TMA-43_scene10', 'JE-TMA-62_scene10','JE-TMA-60_scene10',\n",
    "       'JE-TMA-43_scene11', 'JE-TMA-60_scene11', 'JE-TMA-62_scene11', 'JE-TMA-43_scene13',\n",
    "       'JE-TMA-62_scene12', 'JE-TMA-43_scene14','JE-TMA-60_scene13', 'JE-TMA-62_scene13']      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_order_r = ls_order[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load original\n",
    "'''\n",
    "s_sample = '20210320_SMTBx2-4_JE-TMA-43_60_62'\n",
    "n_neighbors = 30\n",
    "n_pcs = 19\n",
    "n_markers = n_pcs+1\n",
    "resolution = 0.5\n",
    "results_file = f'{s_sample}_{n_neighbors}neighbors_{n_pcs}pcs_{n_markers}markers_leiden{resolution}.h5ad'\n",
    "adata1 = sc.read_h5ad(results_file) \n",
    "\n",
    "d_cluster = {'14': '14: Basal',\n",
    "'5': '5: T cell',\n",
    "'12': '12: T cell',\n",
    "'10': '10: Myoepithelial',\n",
    "'1': '1: Mesenchymal',\n",
    "'16': '16: Prolif.',\n",
    "'15': '15: Vim+ FB (Bx3)',\n",
    "'11': '11: Vim+ FB (Bx4)',\n",
    "'13': '13: Vim+ FB (Bx2)',\n",
    "'7': '7: HER2++',\n",
    "'9': '9: EGFR+ Basal',\n",
    "'3': '3: HER2+',\n",
    "'8': '8: HER2++, Ecad-',\n",
    "'0': '0: ER+ (Bx4)',\n",
    "'2': '2: ER+, PCNA+ ',\n",
    "'4': '4: ER+, EGFR+ (Bx3)',\n",
    "'6': '6: ER+ (Bx2)'}\n",
    "d_clust_names = dict(zip([item[0] for item in d_cluster.items()],[item[1].split(': ')[1] for item in d_cluster.items()]))\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load\n",
    "s_sample = '20210325_SMTBx2-4_JE-TMA-43_60_62'\n",
    "n_neighbors = 31\n",
    "n_pcs = 17\n",
    "n_markers = n_pcs+1\n",
    "resolution = 0.45\n",
    "results_file = f'{s_sample}_{n_neighbors}neighbors_{n_pcs}pcs_{n_markers}markers_leiden{resolution}.h5ad'\n",
    "adata1 = sc.read_h5ad(results_file) \n",
    "print(results_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if resolution == 0.5:\n",
    "    d_cluster = {'14': '14: Basal','12': '12: T cell','16': '16: Prolif.','7': '7: ER+ (Bx2)','13': '13: Luminal (N.Breast)',\n",
    "    '1': '1: ER+ PCNA+ (T47D)','0': '0: ER+ (Bx4)','15': '15: ER+ CK8++ (Bx4)','4': '4: ER+, EGFR+ (Bx3)','18': '18: ER+, EGFR+ (Bx3)',\n",
    "    '17': '17: (Bx3)','10': '10: FB (Bx4)','11': '11: FB (Bx2)','3': '3: CD44+','9': '9: CD44+', '8': '8: EGFR+ Basal',\n",
    "    '5': '5: HER2++','6': '6: HER2+','2': '2: HER2++, Ecad-',}\n",
    "elif resolution == 0.45:\n",
    "    d_cluster = {'15':'15: Basal',\n",
    "                 '12':'12: T cell',\n",
    "                 '16': '16: prolif.',\n",
    "                 '5':'5: ER+, EGFR+ (Bx3)',\n",
    "                '0':'0: ER+ (Bx4)',\n",
    "                '1':'1: ER+, PCNA+',\n",
    "                '7':'7: ER- (Bx2)',\n",
    "                '9':'9: ER+ (Bx2)',\n",
    "                '8':'8: EGFR+ Basal',\n",
    "                '4':'4: HER2+',\n",
    "                '3':'3: HER2+',\n",
    "                '6':'6: HER2+, Ecad-',\n",
    "                '2':'2: Mesenchymal',\n",
    "                '10':'10: Mesenchymal',\n",
    "                '14':'14: fibroblast',\n",
    "                '11':'11: fibroblast',\n",
    "                '13':'13: fibroblast'}\n",
    "d_clust_names = dict(zip([item[0] for item in d_cluster.items()],[item[1].split(': ')[1] for item in d_cluster.items()]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "#sns.set(font_scale=1.19)\n",
    "#seaborn clustermap\n",
    "df_p = pd.DataFrame(data=adata1.raw.X,index=adata1.obs.index,columns=adata1.var.index)\n",
    "df_p['leiden'] = adata1.obs['leiden']\n",
    "g = sns.clustermap(df_p.groupby('leiden').mean().rename({'eccentricity':'eccen.'},axis=1).rename(d_cluster, axis=0),\n",
    "                   z_score=1,figsize=(6.2,6),cmap='viridis',\n",
    "                   vmin=-2,vmax=2,cbar_pos=(.05, .89, .10, .05),cbar_kws={'orientation': 'horizontal','label':'Z-score'}) #(left, bottom, width, height),\n",
    "g.savefig(f'./{s_date}/clustermap_leiden_{resolution}_{n_markers}.pdf',dpi=300)\n",
    "g.savefig(f'./{s_date}/clustermap_leiden_{resolution}_{n_markers}.png',dpi=300)\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": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# stacked bar vertical\n",
    "\n",
    "df = pd.DataFrame(data=adata1.raw.X,index=adata1.obs.index,columns=adata1.var.index)\n",
    "df[f'leiden'] = [int(item) for item in adata1.obs.leiden]\n",
    "s_markers = n_markers\n",
    "k=resolution\n",
    "\n",
    "df['slide'] = [item.split('_')[0] for item in df.index]\n",
    "df['slide_scene'] = [item.split('_cell')[0] for item in df.index]\n",
    "df['slide_scene'] = df.slide_scene.replace({'SMTBx2-5_scene001':'Bx2', 'SMTBx2-5_scene002':'Bx2',\n",
    "       'SMTBx3_scene004':'Bx3', 'SMTBx4-3_scene001':'Bx4',\n",
    "       'SMTBx4-3_scene002':'Bx4'})#.replace(d_order)\n",
    "df['scene'] = [item.split('_')[1] for item in df.index]\n",
    "df_prop = (df.groupby([f'leiden','slide_scene']).PCNA.count())/(df.groupby(['slide_scene']).PCNA.count())\n",
    "df_prop = df_prop.unstack().fillna(value=0).T\n",
    "\n",
    "fig,ax=plt.subplots(figsize=(5,6), dpi=200)\n",
    "df_prop['slide'] =[item.split('_')[0] for item in df_prop.index]\n",
    "#df_prop['scene'] =[item.split('_')[1] for item in df_prop.index]\n",
    "df_prop = df_prop.loc[ls_order_r]\n",
    "df_prop.columns = [str(item) for item in df_prop.columns]\n",
    "#df_prop.rename(d_order).rename(d_cluster,axis=1).plot(kind='barh',stacked=True,ax=ax,legend=True,cmap='tab20',width=0.9)\n",
    "df_prop.plot(kind='barh',stacked=True,ax=ax,legend=True,cmap='tab20',width=0.9)\n",
    "ax.legend(bbox_to_anchor=(1.05, 1.00),ncol=1, fancybox=True,title='Cluster ID')\n",
    "ax.set_xlabel('Fraction of Cells')\n",
    "ax.set_ylabel('Tissue')\n",
    "ax.set_title('')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./{s_date}/StackedBar_{s_markers}markers_{k}Clusters_vertical.pdf')\n",
    "fig.savefig(f'./{s_date}/StackedBar_{s_markers}markers_{k}Clusters_vertical.png')\n",
    "#plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save the cluster ID, not hte annotation\n",
    "#df_prop.to_csv(f'{s_sample}_{n_markers}markers_leiden{resolution}_frac_pos.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.ticker as tic\n",
    "#SMT\n",
    "fig,ax=plt.subplots(figsize=(2.8,3.2),dpi=200)\n",
    "df_plot = df_prop.loc[['Bx2','Bx3','Bx4'],df_prop.dtypes=='float64'].T[::-1]\n",
    "df_plot.plot(kind='barh',ax=ax,legend=True,width=.9)\n",
    "ax.legend(title='Bx', loc='upper left',fancybox=True,borderpad=.2,bbox_to_anchor=(1.05, 1.05))\n",
    "ax.set_xlabel('Fraction of Cells')\n",
    "ax.set_ylabel('')\n",
    "fig.suptitle(f'Cluster Composition: Biopsies',x=.5, y=.92)\n",
    "for tick in ax.yaxis.get_major_ticks():\n",
    "    tick.tick1line.set_markersize(0)\n",
    "    tick.tick2line.set_markersize(0)\n",
    "temp = tic.LinearLocator(numticks=18)\n",
    "ax.yaxis.set_minor_locator(temp)\n",
    "plt.grid(b=True, which='minor', axis='y')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./{s_date}/Barplot_SMT{s_markers}_K{k}.pdf')\n",
    "fig.savefig(f'./{s_date}/Barplot_SMT{s_markers}_K{k}.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_order = ['Bx2', 'Bx3', 'Bx4','AU565-2','AU565-3', 'AU565-4', 'BT474-2','BT474-3', 'BT474-4', \n",
    "       'HCC1143-2', 'HCC1143-4', 'HCC3153-2', 'HCC3153-4', #'JE-TMA-43_scene01','JE-TMA-62_scene01', 'JE-TMA-43_scene10',\n",
    "        'MDAMB-436-2','MDAMB-436-3', 'MDAMB-436-4', 'T47D-2','T47D-3', 'T47D-4',\n",
    "       'N.Breast-2', 'N.Breast-4', 'tonsil-2', 'tonsil-4']\n",
    "d_order = {#'\n",
    "       'JE-TMA-43_scene02':'HCC1143-2', 'JE-TMA-62_scene02':'HCC1143-4',\n",
    "       'JE-TMA-43_scene03':'HCC3153-2', 'JE-TMA-62_scene03':'HCC3153-4', 'JE-TMA-43_scene04':'N.Breast-2',\n",
    "       'JE-TMA-62_scene04':'N.Breast-4', 'JE-TMA-43_scene05':'T47D-2', 'JE-TMA-62_scene05':'T47D-4',\n",
    "       'JE-TMA-43_scene06':'T47D-2', 'JE-TMA-62_scene06':'T47D-4', 'JE-TMA-43_scene07':'tonsil-2',\n",
    "       'JE-TMA-62_scene07':'tonsil-4', 'JE-TMA-43_scene08':'BT474-2', 'JE-TMA-62_scene08':'BT474-4',\n",
    "       'JE-TMA-43_scene09':'BT474-2', 'JE-TMA-62_scene09':'BT474-4',  'JE-TMA-43_scene10':'AU565-2','JE-TMA-62_scene10':'AU565-4',\n",
    "       'JE-TMA-43_scene11':'AU565-2', 'JE-TMA-62_scene11':'AU565-4', 'JE-TMA-43_scene13':'MDAMB-436-2',\n",
    "       'JE-TMA-62_scene12':'MDAMB-436-4', 'JE-TMA-43_scene14':'MDAMB-436-2', 'JE-TMA-62_scene13':'MDAMB-436-4',\n",
    "       'JE-TMA-60_scene13':'MDAMB-436-3', 'JE-TMA-60_scene11':'AU565-3', 'JE-TMA-60_scene10':'AU565-3',\n",
    "       'JE-TMA-60_scene09':'BT474-3', 'JE-TMA-60_scene08':'BT474-3', 'JE-TMA-60_scene06':'T47D-3'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#stacked bar vertical tissue\n",
    "df['coreID'] = df.slide_scene.replace(d_order)\n",
    "df['celltype'] = df.leiden.astype('str').replace(d_clust_names)\n",
    "df_prop = (df.groupby([f'celltype','coreID']).PCNA.count())/(df.groupby(['coreID']).PCNA.count())\n",
    "df_prop = df_prop.unstack().fillna(value=0).T\n",
    "\n",
    "fig,ax=plt.subplots(figsize=(5,3.7), dpi=200)\n",
    "df_prop['slide'] =[item.split('_')[0] for item in df_prop.index]\n",
    "ls_order_r = ls_order[::-1]\n",
    "df_prop = df_prop.loc[ls_order_r]\n",
    "df_prop.columns = [str(item) for item in df_prop.columns]\n",
    "df_prop.plot(kind='barh',stacked=True,ax=ax,legend=True,cmap='tab20',width=0.9) #.rename(d_order).rename(d_clust_names,axis=1)\n",
    "ax.legend(loc='upper left', bbox_to_anchor=(1.1,1.02),ncol=1, fancybox=True,title='Cluster Annotation')\n",
    "ax.set_xlabel('Fraction of Cells')\n",
    "ax.set_ylabel('Tissue')\n",
    "ax.set_title('Cluster Composition: Biopsies Plus Controls')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./{s_date}/StackedBar_{s_markers}markers_{k}Clusters_withcontrols_vert.pdf')\n",
    "fig.savefig(f'./{s_date}/StackedBar_{s_markers}markers_{k}Clusters_withcontrols_vert.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#stacked bar horizontal\n",
    "df['coreID'] = df.slide_scene.replace(d_order)\n",
    "df['celltype'] = df.leiden.astype('str').replace(d_clust_names)\n",
    "df_prop = (df.groupby([f'celltype','coreID']).PCNA.count())/(df.groupby(['coreID']).PCNA.count())\n",
    "df_prop = df_prop.unstack().fillna(value=0).T\n",
    "\n",
    "fig,ax=plt.subplots(figsize=(10,2.5), dpi=200)\n",
    "df_prop['slide'] =[item.split('_')[0] for item in df_prop.index]\n",
    "#df_prop['scene'] =[item.split('_')[1] for item in df_prop.index]\n",
    "df_prop = df_prop.loc[ls_order]\n",
    "df_prop.columns = [str(item) for item in df_prop.columns]\n",
    "df_prop.plot(kind='bar',stacked=True,ax=ax,legend=True,cmap='tab20',width=0.9) #.rename(d_order).rename(d_clust_names,axis=1)\n",
    "ax.legend(loc='upper center', bbox_to_anchor=(1.5, 1.05),ncol=2, fancybox=True,title='Cluster Annotation')\n",
    "ax.set_ylabel('Fraction of Cells')\n",
    "ax.set_xlabel('Tissue')\n",
    "ax.set_title('')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./{s_date}/StackedBar_{s_markers}markers_{k}Clusters_withcontrols.pdf')\n",
    "fig.savefig(f'./{s_date}/StackedBar_{s_markers}markers_{k}Clusters_withcontrols.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot all groups spatially  \n",
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "newcmap = ListedColormap(mpl.cm.tab20.colors)#ListedColormap(mpl.cm.tab20b.colors + mpl.cm.tab20c.colors)\n",
    "from mplex_image import analyze\n",
    "df_pos = analyze.celltype_to_bool(df_p,'leiden')\n",
    "df_xy = df_mi.loc[df_pos.index]\n",
    "ls_scene = ['SMTBx2-5_scene001', 'SMTBx3_scene004', 'SMTBx4-3_scene001', 'SMTBx4-3_scene002']\n",
    "#ls_scene = ['JE-TMA-62_scene04', 'JE-TMA-43_scene04','JE-TMA-62_scene07','JE-TMA-43_scene07']\n",
    "for s_slide in ls_scene:\n",
    "    fig,ax = plt.subplots(figsize=(10,10),dpi=200) #10,10\n",
    "    #plot negative cells\n",
    "    df_scene = df_xy[df_xy.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",
    "        s_color = str(s_color_int)\n",
    "        if len(df_xy[(df_xy.slide_scene==s_slide) & (df_pos.loc[:,s_color])])>=1:\n",
    "            #plot positive cells\n",
    "            ax.scatter(data=df_xy[(df_xy.slide_scene==s_slide) & (df_pos.loc[:,s_color])],x='DAPI_X',y='DAPI_Y',\n",
    "                                                                        label=f'{s_color}',s=0.1,color=newcmap.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) \n",
    "    fig.savefig(f'{codedir}/paper_data/GatingPlots/{s_slide}_clustering_scatterplot.png')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(f'{s_sample}_{n_markers}markers_leiden{resolution}.csv'):\n",
    "    print('saving csv')\n",
    "    df.to_csv(f'{s_sample}_{n_markers}markers_leiden{resolution}.csv')\n",
    "    df_prop.to_csv(f'{s_sample}_{n_markers}markers_leiden{resolution}_frac_pos.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f'{s_sample}_{n_markers}markers_leiden{resolution}.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f'{s_sample}_{n_markers}markers_leiden{resolution}_frac_pos.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.9.5",
   "language": "python",
   "name": "python3.9.5"
  },
  "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
}