DoubleApplication_K157.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Double Antibody Application to a Single Tissue\n",
    "\n",
    "**Question:** How does the stain sensitivity and specificity change when you apply the same antibody panel twice to a single tissue?\n",
    "\n",
    "**Samples:** \n",
    "- A TMA containing HER2+ breast cancer tissue, purchased from Biomax. https://www.biomax.us/tissue-arrays/Breast/Her2B\n",
    "\n",
    "**Method**: We utilized the napari image viewer to overlay all of the markers, set thresholds and create masks. For positive signal, we set a threshold (recorded in our 20211007_napari.py script) that created a mask including all of the positive staining plus any bright artifacts. We then manually erased areas of the mask covering imaging or staining artifacts and saved the mask for future reproducibility. For background signal, we manually selected six regions of the image exhibiting background caused by the artifacts listed above, and again saved the mask. We then extracted the mean intensity of the positive and background areas of the image and used these to calculate SBR.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load libraries\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import sys\n",
    "import numpy as np\n",
    "import os\n",
    "from scipy.stats.stats import pearsonr\n",
    "import itertools\n",
    "from scipy import stats\n",
    "import seaborn as sns\n",
    "\n",
    "matplotlib.rcParams.update({'font.size': 12})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K157 Analysis Steps\n",
    "\n",
    "1 - Estimate dynamic range and first/second correlation in ROIs using percentile \n",
    "\n",
    "2 - Visualize intensity distribution of early (i.e. first application) and late (i.e. second application) signal and background\n",
    "\n",
    "3 - Use napari to manually select foreground/ background pixels. Masks are saved for reproducibility. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set location of files\n",
    "#os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cycIF_ValidationStudies/cycIF_Validation')\n",
    "codedir = os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load positive (thresholded using AF subtracted)\n",
    "df_pos = pd.read_csv(f'{codedir}/Data/Her2B-K157-6S_DAPI22Q_Nuclei1000_ManualPositive.txt',index_col=0,sep='\\t')\n",
    "df_pos['slide'] = [item.split('_')[0] for item in df_pos.index]\n",
    "df_pos['scene'] = [item.split('_')[1] for item in df_pos.index]\n",
    "\n",
    "#OLD \n",
    "'''read in full features (all locations, cells) (too big)\n",
    "#df_mi = pd.read_csv(f'{codedir}/BigData/features_Her2B-K157-6S_MeanIntensity.tsv',index_col=0,sep='\\t')\n",
    "\n",
    "#df_mi['scene'] = [item.split('_')[0] for item in df_mi.index]# save smaller dataframes\n",
    "#for s_scene in sorted(set(df_mi.scene)):\n",
    "#    df_mi[df_mi.scene==s_scene].to_csv(f'{codedir}/Data/features_Her2B-K157-6S-{s_scene}_MeanIntensity.csv')\n",
    "#still too big\n",
    "#s_scene = 'scene002'\n",
    "#df_mi.loc[df_mi.scene==s_scene,~df_mi.columns.str.contains('_Cells')].to_csv(f'{codedir}/Data/features_Her2B-K157-6S-{s_scene}_MeanIntensity.csv')\n",
    "'''\n",
    "#load \n",
    "df_mi = pd.DataFrame()\n",
    "for s_scene in ['scene002', 'scene004', 'scene005', 'scene006', 'scene008', 'scene009']:\n",
    "    df_mi = df_mi.append(pd.read_csv(f'{codedir}/Data/features_Her2B-K157-6S-{s_scene}_MeanIntensity.csv',index_col=0))\n",
    "#Define tissue compartments\n",
    "\n",
    "#define tumor and immune cell types\n",
    "ls_tumor = ['CK19_Ring','CK5_Ring','Ecad_Ring','CK14_Ring'] # #'CK7_Ring'= poor quench, CK17_Ring = too many pos comapred to CK5 and CK14\n",
    "ls_exclude = ['CD31_Ring','CD68_Ring','CD4_Ring'] \n",
    "df_pos['all-tumor'] = df_pos.loc[:,ls_tumor].sum(axis=1)>=1 \n",
    "df_pos['not-tumor'] = df_pos.loc[:,ls_exclude].sum(axis=1)>=1\n",
    "df_pos['not-tumor-not-Tcell'] = df_pos.loc[:,['CD31_Ring','CD68_Ring',]].sum(axis=1)>=1\n",
    "df_pos['tumor'] = (df_pos.loc[:,'all-tumor'].astype(int) -  df_pos.loc[:,'not-tumor'].astype(int))==1\n",
    "df_pos['Endothelial'] = df_pos.loc[:,'CD31_Ring']\n",
    "\n",
    "df_pos['Macrophage'] = (df_pos.loc[:,'CD68_Ring'].astype(int) - df_pos.loc[:,'CD31_Ring'].astype(int)) == 1\n",
    "\n",
    "\n",
    "df_pos['Tcell'] = (df_pos.loc[:,'CD4_Ring'].astype(int) -  df_pos.loc[:,'not-tumor-not-Tcell'].astype(int))==1\n",
    "df_pos['Stromal'] = (df_pos.loc[:,['tumor','Tcell','Macrophage','Endothelial']].sum(axis=1) ==0)\n",
    "df_pos = df_pos.drop(['all-tumor','not-tumor','not-tumor-not-Tcell'],axis=1)\n",
    "\n",
    "#check\n",
    "#df_pos.loc[:,['tumor','Tcell','Mac','Endothelial','Stromal']].sum().sum()\n",
    "#len(df_pos)\n",
    "#(df_pos.loc[:,['tumor','Tcell','Mac','Endothelial','Stromal']].sum(axis=1)!=1).sum()\n",
    "#five compartments\n",
    "df_pos.loc[:,['tumor','Tcell','Macrophage','Endothelial','Stromal']].sum()\n",
    "#Annotate the compartment of signal and background for each marker\n",
    "#dd key is forground, d key is cell type, d item is bg\n",
    "dd_exclusive = {'Tcell':{ 'FoxP3_Nuclei':'tumor',\n",
    "    'GRNZB_Nuclei':'tumor',\n",
    "    'PD1_Ring':'tumor',\n",
    "    'CD20_Ring':'tumor',\n",
    "    'CD45_Ring':'tumor',\n",
    "    'CD8_Ring':'tumor',\n",
    "    'CD4_Ring':'tumor'},\n",
    "                \n",
    "'tumor':{'CK8_Ring':'Stromal',\n",
    "    'CK5_Ring':'Stromal',\n",
    "    'CK17_Ring':'Stromal',\n",
    "    'HER2_Ring':'Stromal',\n",
    "    'Ecad_Ring':'Stromal',\n",
    "    'CK14_Ring':'Stromal',\n",
    "    'ER_Nuclei':'Stromal',\n",
    "    'PgR_Nuclei':'Stromal',\n",
    "    'CK19_Ring':'Stromal'},\n",
    "'Endothelial': {'CD31_Ring':'tumor',},\n",
    "'Macrophage':{'CD68_Ring':'tumor',},\n",
    "'Stromal':{'ColIV_Ring':'tumor',\n",
    "    'aSMA_Ring':'Tcell',\n",
    "    'PDPN_Ring':'tumor',\n",
    "    #'ColI_Ring':'tumor', #different lot of antibody between early and late\n",
    "    'CD44_Ring':'Endothelial',\n",
    "    'Vim_Ring':'tumor',}\n",
    "               }\n",
    "d_scene = {'AR_Nuclei':['scene008'],\n",
    "           'aSMA_Ring':['scene002','scene006'],\n",
    "           'CD4_Ring':['scene002','scene008'],\n",
    "           'CD8_Ring':['scene002','scene008'],\n",
    "           'CD20_Ring':['scene008'], #'scene002',\n",
    "           'CD31_Ring':['scene002','scene006','scene008'],\n",
    "           'CD44_Ring':['scene002','scene006','scene008'],\n",
    "           'CD45_Ring':['scene002','scene008'],\n",
    "           'CD68_Ring':['scene002','scene006','scene008'],\n",
    "           'CK5_Ring':['scene006'],\n",
    "           'CK7_Ring':['scene002','scene006'],\n",
    "           'CK8_Ring':['scene002','scene006','scene008'],\n",
    "           'CK14_Ring':['scene006'],\n",
    "           'CK17_Ring':['scene006'],\n",
    "           'CK19_Ring':['scene002','scene006','scene008'],\n",
    "           'ColI_Ring':['scene006'],\n",
    "           'ColIV_Ring':['scene006'],\n",
    "           'cPARP_Nuclei':['scene002'],\n",
    "           'Ecad_Ring':['scene002','scene008'],\n",
    "           'ER_Nuclei':['scene008'],\n",
    "           'FoxP3_Nuclei':['scene002','scene008'],\n",
    "           'GRNZB_Nuclei':['scene002','scene006'],\n",
    "           'H3K4_Nuclei':['scene002','scene008'],\n",
    "           'H3K27_Nuclei':['scene002','scene008'],\n",
    "           'HER2_Ring':['scene002','scene006'],\n",
    "           'Ki67_Nuclei':['scene002','scene006','scene008'],\n",
    "           'LamAC_Nuclei':['scene002','scene006','scene008'],\n",
    "           'LamB1_Nuclei':['scene002','scene008'],\n",
    "           'LamB2_Nuclei':['scene002','scene006','scene008'],\n",
    "           'PCNA_Nuclei':['scene002','scene008'],\n",
    "           'PD1_Ring':['scene008'],\n",
    "           'PDPN_Ring':['scene002','scene008'],\n",
    "           'PgR_Nuclei':['scene008'],\n",
    "           'pRB_Nuclei':['scene008'],\n",
    "           'pS6RP_Ring':['scene008'],\n",
    "           'S100_Nuclei':['scene006'],\n",
    "           'Vim_Ring':['scene002','scene006','scene008'],\n",
    "          }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# dynamic range, SBR\n",
    "\n",
    "compare first application to second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ROIs that have known positive cells\n",
    "d_crop ={'Her2B-K157_scene002':(2160,4350),\n",
    "    'Her2B-K157_scene006':(4500,4100),\n",
    "    'Her2B-K157_scene008':(2800,950),}\n",
    "\n",
    "\n",
    "df_xy = df_pos.loc[:,['DAPI_X','DAPI_Y','scene','slide']]\n",
    "df_xy['slide'] = [item.split(\"-6S\")[0] for item in df_xy.slide] \n",
    "df_xy['slide_scene'] = df_xy.slide + '_' + df_xy.scene\n",
    "\n",
    "df_rois = pd.DataFrame()\n",
    "for s_scene, t_coord in d_crop.items():\n",
    "        x_min = t_coord[0]\n",
    "        x_max = t_coord[0] + 800\n",
    "        y_min = t_coord[1]\n",
    "        y_max = t_coord[1] + 800\n",
    "        df_scene = df_xy[df_xy.slide_scene==s_scene]\n",
    "        df_scene = df_scene.loc[(x_max > df_scene.DAPI_X) & (df_scene.DAPI_X > x_min)]\n",
    "        df_scene = df_scene.loc[(y_max > df_scene.DAPI_Y) & (df_scene.DAPI_Y > y_min)]\n",
    "        df_rois = df_rois.append(df_mi.loc[[item.split(\"6S_\")[1] for item in df_scene.index],:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dynamic range\n",
    "upper_quantile = 0.995\n",
    "#calculate 5th and 98th percentile to estimate dynamic range\n",
    "df_result = pd.DataFrame()\n",
    "df_early = pd.DataFrame()\n",
    "df_late = pd.DataFrame()\n",
    "\n",
    "for s_marker, ls_scene in d_scene.items():   \n",
    "        for s_scene in ls_scene:\n",
    "            se_early = df_rois.loc[df_mi.scene==s_scene,s_marker]\n",
    "            se_late = df_rois.loc[df_mi.scene==s_scene,f'{s_marker.split(\"_\")[0]}b_{s_marker.split(\"_\")[1]}']\n",
    "            #early\n",
    "            s_index = f'{s_scene}_{s_marker}'\n",
    "            df_early.loc[s_index,'low05'] = se_early.quantile(0.04)\n",
    "            df_early.loc[s_index,'high95'] = se_early.quantile(upper_quantile)\n",
    "            df_early.loc[s_index,'order'] = 'early'\n",
    "            df_early.loc[s_index,'core'] = s_scene\n",
    "            df_early.loc[s_index,'marker'] = s_marker.split('_')[0]\n",
    "            #late\n",
    "            df_late.loc[s_index,'low05'] = se_late.quantile(0.04)\n",
    "            df_late.loc[s_index,'high95'] = se_late.quantile(upper_quantile)\n",
    "            df_late.loc[s_index,'order'] = 'late'\n",
    "            df_late.loc[s_index,'core'] = s_scene\n",
    "            df_late.loc[s_index,'marker'] = s_marker.split('_')[0]\n",
    "            # coorelation\n",
    "            r, p  = stats.pearsonr(se_early, se_late)\n",
    "            df_result.loc[s_index,'pearson'] = r\n",
    "            df_result.loc[s_index,'pval'] = p\n",
    "            df_result.loc[s_index,'core'] = s_scene\n",
    "            df_result.loc[s_index,'marker'] = s_marker.split('_')[0]\n",
    "            \n",
    "df_early['dynamic_range'] = df_early.high95 - df_early.low05\n",
    "df_early['estimated_SBR'] = (df_early.high95- df_early.low05)/df_early.low05\n",
    "\n",
    "df_late['dynamic_range'] = df_late.high95 - df_late.low05\n",
    "df_late['estimated_SBR'] = (df_late.high95- df_late.low05)/df_late.low05\n",
    "\n",
    "ls_col = df_early.columns[df_early.dtypes=='float64']\n",
    "df_ratio = df_early.loc[:,ls_col]/df_late.loc[:,ls_col]\n",
    "\n",
    "df_result = df_result.merge(df_ratio, left_index=True, right_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dynamic range\n",
    "df_order = pd.DataFrame()\n",
    "for s_marker in df_result.groupby('marker').dynamic_range.mean().sort_values().index:\n",
    "    df_order = df_order.append(df_result[df_result.marker==s_marker])\n",
    "fig, ax=plt.subplots(figsize=(10,4),dpi=200)\n",
    "sns.lineplot(\n",
    "    data=df_order, x=\"marker\", y=\"dynamic_range\", err_style=\"bars\", ci='sd', ax=ax\n",
    ")\n",
    "ax.axhline(1,color='black')\n",
    "labels=ax.get_xticklabels()\n",
    "ax.tick_params(axis='x', labelrotation= 90) \n",
    "ax.set_ylabel(f'Dynamic Range 0.04 - {upper_quantile}')\n",
    "ax.set_title('Ratio of First versus Second Application: Dynamic Range')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/K157_dynamic_range{upper_quantile}_ratio.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#percent antibodies earlier higher dynmic range\n",
    "sum(df_order.groupby('marker').dynamic_range.mean() > 1)/len(df_order.groupby('marker').dynamic_range.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#percent antibodies less than 50% change\n",
    "sum((df_order.groupby('marker').dynamic_range.mean() < 1.5) & (df_order.groupby('marker').dynamic_range.mean() > .5))/len(df_order.groupby('marker').dynamic_range.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# correlation\n",
    "df_order = pd.DataFrame()\n",
    "for s_marker in df_result.groupby('marker').pearson.mean().sort_values().index:\n",
    "    df_order = df_order.append(df_result[df_result.marker==s_marker])\n",
    "fig, ax=plt.subplots(figsize=(10,4),dpi=200)\n",
    "sns.lineplot(\n",
    "    data=df_order, x=\"marker\", y=\"pearson\", err_style=\"bars\", ci='sd', ax=ax\n",
    ")\n",
    "ax.axhline(1,color='black')\n",
    "labels=ax.get_xticklabels()\n",
    "ax.set_ylabel(\"Pearson's r\")\n",
    "ax.tick_params(axis='x', labelrotation= 90) \n",
    "ax.set_title('First versus Second Application: Pearson Correlation')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/K157_pearson.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#percent antibodies high pearson\n",
    "sum(df_order.groupby('marker').pearson.mean() > .8)/len(df_order.groupby('marker').pearson.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SBR (use only scenes that are positive)\n",
    "df_order = pd.DataFrame()\n",
    "for s_marker in df_result.groupby('marker').estimated_SBR.mean().sort_values().index:\n",
    "    df_order = df_order.append(df_result[df_result.marker==s_marker])\n",
    "fig, ax=plt.subplots(figsize=(10,4),dpi=200)\n",
    "sns.lineplot(\n",
    "    data=df_order, x=\"marker\", y=\"estimated_SBR\", err_style=\"bars\", ci='sd', ax=ax\n",
    ")\n",
    "ax.axhline(1,color='black')\n",
    "labels=ax.get_xticklabels()\n",
    "ax.set_ylabel(\"Est. SBR\")\n",
    "ax.tick_params(axis='x', labelrotation= 90) \n",
    "ax.set_title('Ratio of First versus Second Application: Estimated SBR')\n",
    "ax.set_ylim(0,15)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/K157_SBR_ratio.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#percent est SBR earlier higher\n",
    "sum(df_order.groupby('marker').estimated_SBR.mean() > 1)/len(df_order.groupby('marker').estimated_SBR.mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SBR\n",
    "df_order = pd.DataFrame()\n",
    "for s_marker in df_late.groupby('marker').estimated_SBR.mean().sort_values().index:\n",
    "    df_order = df_order.append(df_late[df_late.marker==s_marker])\n",
    "fig, ax=plt.subplots(figsize=(10,4),dpi=200)\n",
    "sns.lineplot(\n",
    "    data=df_order, x=\"marker\", y=\"estimated_SBR\", label='late', ax=ax, err_style=\"bars\", ci='sd',\n",
    ")\n",
    "ax.axhline(1.5,color='black')\n",
    "labels=ax.get_xticklabels()\n",
    "ax.tick_params(axis='x', labelrotation= 90) \n",
    "ax.set_title('Second Application SBR')\n",
    "ax.get_legend().remove()\n",
    "ax.set_ylim(0,15)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/K157_lateSBR.png')\n",
    "\n",
    "df_order = pd.DataFrame()\n",
    "for s_marker in df_early.groupby('marker').estimated_SBR.mean().sort_values().index:\n",
    "    df_order = df_order.append(df_early[df_early.marker==s_marker])\n",
    "fig, ax=plt.subplots(figsize=(10,4),dpi=200)\n",
    "sns.lineplot(\n",
    "    data=df_order, x=\"marker\", y=\"estimated_SBR\", label='early', ax=ax, err_style=\"bars\", ci='sd',\n",
    ")\n",
    "ax.axhline(1.5,color='black')\n",
    "labels=ax.get_xticklabels()\n",
    "ax.tick_params(axis='x', labelrotation= 90) \n",
    "ax.set_title('First Application SBR')\n",
    "ax.get_legend().remove()\n",
    "ax.set_ylim(0,35)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/K157_earlySBR.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sum(df_early.groupby('marker').estimated_SBR.mean() > 1.5)/len(df_early.groupby('marker').estimated_SBR.mean()))\n",
    "print(sum(df_late.groupby('marker').estimated_SBR.mean() > 1.5)/len(df_late.groupby('marker').estimated_SBR.mean()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#produce single cell data/histograms\n",
    "\n",
    "for s_fg, d_exclusive in dd_exclusive.items():\n",
    "    fg_index = [item.split('Her2B-K157-6S_')[1] for item in  df_pos[df_pos.loc[:,s_fg]].index]\n",
    "    for s_marker, s_bg in d_exclusive.items():    \n",
    "        fig, ax = plt.subplots(2,1, figsize=(3,3.5),sharex=True,dpi=200)\n",
    "        ax=ax.ravel()\n",
    "        #early\n",
    "        bg_index = [item.split('Her2B-K157-6S_')[1] for item in  df_pos[df_pos.loc[:,s_bg]].index]\n",
    "        se_fg = df_mi.loc[fg_index,s_marker]\n",
    "        se_bg = df_mi.loc[bg_index,f\"{s_marker.split('_')[0]}_{s_marker.split('_')[1].replace('Ring','Nucleus').replace('Nuclei','Ring').replace('Nucleus','Nuclei')}\"]\n",
    "        sns.distplot(se_fg,bins=500,ax=ax[0],kde=False,label='S',hist_kws=dict(alpha=0.4))                       \n",
    "        sns.distplot(se_bg,bins=500,ax=ax[0],kde=False,label='BG',hist_kws=dict(alpha=0.4))\n",
    "        ax[0].legend()\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}')\n",
    "        ax[0].set_xlabel('Early')\n",
    "        #late\n",
    "        se_fg = df_mi.loc[fg_index,f'{s_marker.split(\"_\")[0]}b_{s_marker.split(\"_\")[1]}']\n",
    "        se_bg = df_mi.loc[bg_index,f\"{s_marker.split('_')[0]}b_{s_marker.split('_')[1].replace('Ring','Nucleus').replace('Nuclei','Ring').replace('Nucleus','Nuclei')}\"]\n",
    "        sns.distplot(se_fg,bins=500,ax=ax[1],kde=False,label='S',hist_kws=dict(alpha=0.4))                       \n",
    "        sns.distplot(se_bg,bins=500,ax=ax[1],kde=False,label='BG',hist_kws=dict(alpha=0.4))\n",
    "        ax[1].legend()\n",
    "        ax[1].set_xlabel('Late')\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/Her2B-K157/K157_histogram_{s_marker.split(\"_\")[0]}.png')\n",
    "        if s_marker != 'CD45_Ring':\n",
    "            plt.close(fig)\n",
    "#kde \n",
    "\n",
    "for s_fg, d_exclusive in dd_exclusive.items():\n",
    "    fg_index = [item.split('Her2B-K157-6S_')[1] for item in  df_pos[df_pos.loc[:,s_fg]].index]\n",
    "    for s_marker, s_bg in d_exclusive.items():    \n",
    "        fig, ax = plt.subplots(2,1, figsize=(3,3.5),dpi=200)\n",
    "        ax=ax.ravel()\n",
    "        #early\n",
    "        bg_index = [item.split('Her2B-K157-6S_')[1] for item in  df_pos[df_pos.loc[:,s_bg]].index]\n",
    "        se_fg = df_mi.loc[fg_index,s_marker]\n",
    "        se_bg = df_mi.loc[bg_index,f\"{s_marker.split('_')[0]}_{s_marker.split('_')[1].replace('Ring','Nucleus').replace('Nuclei','Ring').replace('Nucleus','Nuclei')}\"]\n",
    "        i_max = np.quantile(se_fg.append(se_bg),.998)\n",
    "        sns.kdeplot(se_fg,ax=ax[0],shade=True,label='S',alpha=0.4,clip=(0,i_max))                       \n",
    "        sns.kdeplot(se_bg,ax=ax[0],shade=True,label='BG',alpha=0.4,clip=(0,i_max))\n",
    "        ax[0].legend()\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}')\n",
    "        ax[0].set_yticklabels('')\n",
    "        ax[0].set_xlabel('Early')\n",
    "        ax[0].set_xlim(0,i_max)\n",
    "        #late\n",
    "        se_fg = df_mi.loc[fg_index,f'{s_marker.split(\"_\")[0]}b_{s_marker.split(\"_\")[1]}']\n",
    "        se_bg = df_mi.loc[bg_index,f\"{s_marker.split('_')[0]}b_{s_marker.split('_')[1].replace('Ring','Nucleus').replace('Nuclei','Ring').replace('Nucleus','Nuclei')}\"]\n",
    "        sns.kdeplot(se_fg,ax=ax[1],shade=True,label='S',alpha=0.4,clip=(0,i_max))                       \n",
    "        sns.kdeplot(se_bg,ax=ax[1],shade=True,label='BG',alpha=0.4,clip=(0,i_max))\n",
    "        ax[1].legend()\n",
    "        ax[1].set_yticklabels('')\n",
    "        ax[1].set_xlabel('Late')\n",
    "        ax[1].set_xlim(0,i_max)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/Her2B-K157/K157_kde_{s_marker.split(\"_\")[0]}.png')\n",
    "        if s_marker != 'CD45_Ring':\n",
    "            plt.close(fig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Napari Analysis\n",
    "\n",
    "Using the interactive image viewer for Python, [napari](https://napari.org/), we manually selected forground and background pixels to quantify SBR in a targeted method. Napari code runs on in-house generated multi-page ometiffs. For each ROI, the foreground and background pixel intensity was saved to csv, which are further analyzed here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lead fg and bg data\n",
    "df_data = pd.DataFrame()\n",
    "s_sample = 'Her2B-K157'\n",
    "for s_dir in os.listdir('./Data'):\n",
    "    if s_dir.find(f'result_{s_sample}') > -1:\n",
    "        print(s_dir)\n",
    "        df_data.loc[s_dir,'slide'] = s_sample\n",
    "        df_data['scene'] = [item.split('_')[1] for item in df_data.index]\n",
    "df_data['rounds'] = [item.split('_')[2] for item in df_data.index]\n",
    "df_data['marker'] = [item.split('_')[3] for item in df_data.index]\n",
    "df_data['region'] = [item.split('_')[4].split('.')[0] for item in df_data.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_result = pd.DataFrame()\n",
    "for s_index in df_data.index:\n",
    "    df = pd.read_csv(f'Data/{s_index}',index_col=0)\n",
    "    df['rounds'] = df_data.loc[s_index,'rounds'] \n",
    "    df['marker'] = df_data.loc[s_index,'marker'] \n",
    "    df['region'] = df_data.loc[s_index,'region']\n",
    "    df_result = df_result.append(df)\n",
    "df_result['mean_intensity'] = df_result.mean_intensity*256\n",
    "df_result.reset_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_marker in ['CD8','CD45']:\n",
    "    fig, ax = plt.subplots(figsize=(3,2.5),dpi=150)\n",
    "    sns.boxplot(x='rounds', y='mean_intensity', hue='region', data=df_result[df_result.marker==s_marker], fliersize=0, ax=ax,hue_order=['fg','bg'])\n",
    "    #sns.stripplot(x='panel', y='mean_intensity', hue='region', data=df_result,ax=ax,palette='dark')\n",
    "    ax.legend(bbox_to_anchor=(.85,.35)) #loc=4\n",
    "    ax.set_title(f'{s_marker} Foreground vs. Background')\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/cmIF_{s_sample}_{s_marker}_FGvBG.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate sbr\n",
    "df_mean = df_result.groupby(['rounds','marker','region']).mean()\n",
    "df_bg = df_result.loc[df_result.region=='bg'].copy()\n",
    "for s_index in df_bg.index:\n",
    "    df_bg.loc[s_index,'fg'] =df_mean.loc[(df_bg.loc[s_index,'rounds'],df_bg.loc[s_index,'marker'],'fg'),'mean_intensity']\n",
    "    df_bg['SBR'] = df_bg.fg/df_bg.mean_intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_bg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "for s_marker in ['CD8','CD45']:\n",
    "    fig, ax = plt.subplots(figsize=(3,2.5),dpi=150)\n",
    "    sns.boxplot(x='rounds', y='SBR', data=df_bg[df_bg.marker==s_marker], fliersize=0, ax=ax)\n",
    "    sns.stripplot(x='rounds', y='SBR', data=df_bg[df_bg.marker==s_marker], palette='dark', ax=ax)\n",
    "    x = df_bg[(df_bg.marker==s_marker) & (df_bg.rounds==df_bg[df_bg.marker==s_marker].groupby('rounds').mean().index[0])].SBR\n",
    "    y = df_bg[(df_bg.marker==s_marker) & (df_bg.rounds==df_bg[df_bg.marker==s_marker].groupby('rounds').mean().index[1])].SBR\n",
    "    statistic, pvalue = stats.ttest_ind(x,y)\n",
    "    \n",
    "    ax.set_title(f'{s_sample} {s_marker} \\n p = {pvalue:.2}')\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/cmIF_{s_sample}_{s_marker}_SBR.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": 2
}