SinglevsCyclic_44290.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of Standard IF to Cyclic IF <a name=\"contents\"></a>\n",
    "\n",
    "**Question:** How does the signal-to-background ratio and staining specificity compare in standard IF vs. cyclic IF?\n",
    "\n",
    "**Samples:**\n",
    "\n",
    "[Single vs Cyclic](#svc)\n",
    "- Tissue ID 44290: HER2+/ER+ breast cancer. Section 112 stained with a 5 round cyclic IF protocol, sections 113 to 116 stained with the same antibodies in a standard IF protocol.\n",
    "- Tissue ID 44294: Adjacent normal breast from the patient above. Section 116 is cyclic, section 117 to 120 are standard IF.\n",
    "\n",
    "**Method**: For each stain, pixel intensity was manually thresholded to separate the positive pixels from the negative, in Extended1_single_vs_cyclic.ipynb notebook. We also automatically estimated dynamic range using the 5th and 98th percentile of single cell intensity in the tissue, plotted dyamic range, and compared SBR and estimated dynamic range.\n",
    "\n",
    "[Tissue Loss](#tissue)\n",
    "- Biomax TMA 808L2: multiple tissue tumor and normal, for tissue loss analysis\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "matplotlib.rcParams.update({'font.size': 16})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notes 44290 (HER2+ Breast Cancer)\n",
    "The single stains have much worse edge effect.\n",
    "Define 2 ROIs with similar composition, avoiding edge regions. (x,y,w,h)\n",
    "\n",
    "44290-112: (cyclic)\n",
    "- (3200,5500,1500,1500) ROI1 large tumor nest upper left\n",
    "- (5000,8800,1500,1800) ROI2 large tumor nests lower right\n",
    "- only two regions in 44290 are suitable, due to autofluorescent regions that are different between single and cyclic stains. \n",
    "\n",
    "44290-113 (single PCNA, HER2, ER, CD45)\n",
    "- (800,3050,1500,1500)\n",
    "- (2650,6400,1500,1800)\n",
    "\n",
    "44290-114 (single pHH3, CK14, CD44, CK5)\n",
    "- (900,2900,1500,1500)\n",
    "- (2650,6250,1500,1800)\n",
    "\n",
    "44290-115 (single Vim, CK7, PD1, Lamin A/C)\n",
    "- (900,2950,1500,1500)\n",
    "- (2850,6200,1500,1800)\n",
    "\n",
    "44290-116 (single aSMA, CD68, Ki67, Ecad)\n",
    "- (850,3000,1500,1500)\n",
    "- (2650,6250,1500,1800)\n",
    "\n",
    "## Notes 44294 (Normal Breast)\n",
    "The single stains on normal tissue don't have such bad edge effect, although they do have to autoflourescence in the DAPI negative areas (i.e. middle of the duct). We can use the whole segmented images, as segmentation based on DAPI will ignore most of the worst autofluorescence areas.\n",
    "\n",
    "44294-116 was not cropped before segmentation, while the other tissues were.\n",
    "Therefore, I will use an ROI to make 44294-116 the same size as the other tissues.\n",
    "\n",
    "44294-116:\n",
    "- (7200,3580,12725,12365) (scene001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Data <a name=\"svc\"></a>\n",
    "\n",
    "For single versus cyclic\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "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": [
    "#import dataframes with mean intensity/centroid from segmentation and feature extraction\n",
    "ls_slide =  ['44290-112', '44290-113', '44290-114', '44290-115', '44290-116', \n",
    "             '44294-116', '44294-117', '44294-118', '44294-119', '44294-120']\n",
    "df_mi = pd.DataFrame()\n",
    "for s_slide in ls_slide:\n",
    "    df_mi = df_mi.append(pd.read_csv(f'{codedir}/Data/features_{s_slide}_singlevscyclic.csv',index_col=0))\n",
    "#add the slide/scene information \n",
    "df_mi['scene'] = [item.split('_')[1] for item in df_mi.index]\n",
    "df_mi['slide'] = [item.split('_')[0] for item in df_mi.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load threshold, SNR data\n",
    "rootdir = '/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cycIF_ValidationStudies/cycIF_Validation'\n",
    "df_t = pd.read_csv(f'{rootdir}/Metadata/44290/SNR_single_vs_cyclic.csv',index_col=0)\n",
    "\n",
    "#load metadata from feature extraction\n",
    "df_m = pd.read_csv(f'{rootdir}/Metadata/44290/metadata_single_vs_cyclic.csv',index_col=0)\n",
    "k=1\n",
    "df_sp = pd.DataFrame()\n",
    "for s_thresh in ['minimum','minimum90','minimum110']: #\n",
    "    filename = f'{rootdir}/Data/single_vs_cyclic_SBR_SP{k}_{s_thresh}.csv'\n",
    "    if os.path.exists(filename):\n",
    "        df = pd.read_csv(filename,index_col='filename')\n",
    "        df['threshold_level'] = s_thresh\n",
    "        df_sp = df_sp.append(df)\n",
    "#alternative regions\n",
    "df_alt = pd.DataFrame()\n",
    "for s_type in ['tum','norm']:\n",
    "    df = pd.read_csv(f'{rootdir}/Data/single_vs_cyclic_SBR_SP1_minimum_alt_{s_type}.csv', index_col='filename')\n",
    "    df['tum_norm'] = s_type\n",
    "    df_alt = df_alt.append(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import dataframes with mean intensity/centroid from segmentation and feature extraction\n",
    "#df_mi =pd.read_csv(f'{codedir}/BigData/features_singlevscyclic.csv',index_col=0)\n",
    "#df_mi['slide'] = [item.split('_')[0] for item in df_mi.index]\n",
    "#for s_slide in sorted(set(df_mi.slide)):\n",
    "#    df_mi[df_mi.slide==s_slide].to_csv(f'{codedir}/Data/features_{s_slide}_singlevscyclic.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#annotation\n",
    "\n",
    "df_t['experiment'] = [item.split('_')[0] for item in df_t.condition]\n",
    "df_t['experimenttype'] = [item.split('_')[2] for item in df_t.condition]\n",
    "df_t = df_t.sort_values(['color','rounds'])\n",
    "\n",
    "#calculate actual signal (ie measure signal minus background)\n",
    "df_t['meanpos'] = df_t.meanpos - df_t.meanneg\n",
    "\n",
    "#add info\n",
    "print(len(df_sp))\n",
    "print(len(df_t))\n",
    "df_sp = df_sp.merge(df_t.loc[:,['experiment', 'experimenttype','exposure', 'refexp','minimum']],left_index=True,right_index=True)\n",
    "df_sp = df_sp.reset_index().rename({'level_0':'filename'},axis=1).drop(['Unnamed: 0','index'],axis=1)\n",
    "print(len(df_sp))\n",
    "df_alt = df_alt.merge(df_t.loc[:,['experiment', 'experimenttype','exposure', 'refexp','minimum']],left_index=True,right_index=True)\n",
    "df_alt = df_alt.reset_index().rename({'level_0':'filename'},axis=1).drop(['Unnamed: 0','index'],axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## super pixel SBR\n",
    "\n",
    "\n",
    "the authors talk about “quantification of the signal” using a threshold, but this is not quantification of the signal, it’s a count of cells considered as “positive” without taking in consideration the full range of value of the signal\n",
    "\n",
    "- compare SBR and estimated dynamic range\n",
    "\n",
    "how generalizable thresholds are to other regions of the images\n",
    "\n",
    "- use superpixels to split image into regions\n",
    "\n",
    "- test another ROI in normal breast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate SBR\n",
    "for s_measure in ['mean_intensity','quartiles-0','quartiles-1','quartiles-2','quartiles-3','quartiles-4']:\n",
    "    df_sp.loc[:,f'{s_measure}_SBR'] = df_sp.loc[:,f'{s_measure}_fg']/df_sp.loc[:,'mean_intensity_bg']\n",
    "df_sp['dynamic_range_SBR'] = df_sp.loc[:,'quartiles-4_SBR'] - df_sp.loc[:,'quartiles-0_SBR']\n",
    "df_sp['dynamic_range'] = df_sp.loc[:,'quartiles-4_fg'] - df_sp.loc[:,'quartiles-0_fg'] #\n",
    "df_sp['tissue_label'] = df_sp.experimenttype + '_' + df_sp.label.astype('str')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_sp.columns\n",
    "ls_order_sp = ['LaminAC','CK14','ER','CD45','aSMA','PD1','CD68',\n",
    "               'PCNA', 'Her2','CK5', 'CD44','Vim','Ki67',  'Ecad',\n",
    "             'CK7',  \n",
    "            ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#plot\n",
    "d_measure = {'mean_intensity_SBR':'SBR',\n",
    "    'dynamic_range':'Dynamic Range',\n",
    "    'mean_intensity_fg':'Mean Intensity',\n",
    "            'mean_intensity_bg':'Background Intensity'\n",
    "            }\n",
    "d_result = {}\n",
    "d_plot = {}\n",
    "for tu_level in [('minimum','minimum')]:#\n",
    "    for s_measure, s_name in d_measure.items():   \n",
    "        df_result = pd.DataFrame()\n",
    "        df_plot = pd.DataFrame()\n",
    "        for s_marker in df_sp.marker.unique():\n",
    "            #s_marker = 'CK7'\n",
    "            df_marker = df_sp[df_sp.marker==s_marker].sort_values(by='experimenttype')\n",
    "            if not len(df_marker==2):\n",
    "                print(s_marker)\n",
    "            df = df_marker[df_marker.experiment=='cyclic'].merge(df_marker[df_marker.experiment=='single'],on='tissue_label',suffixes=('_cyclic','_single'))\n",
    "            d_rename = dict(zip([i for i in range(len(df))],[item.split('_')[0] for item in df.tissue_label]))\n",
    "            se_sing = df.loc[df.loc[:,'threshold_level_single']==tu_level[0],f'{s_measure}_single'].rename(d_rename)\n",
    "            se_cyc = df.loc[df.loc[:,'threshold_level_cyclic']==tu_level[1],f'{s_measure}_cyclic'].rename(d_rename)\n",
    "\n",
    "            df_plot[s_marker] = (se_sing/se_cyc)\n",
    "            se_sing.name = f'{s_marker}_single'\n",
    "            se_cyc.name = f'{s_marker}_cyclic'\n",
    "            df_result=df_result.append(se_sing)\n",
    "            df_result=df_result.append(se_cyc)\n",
    "            #break\n",
    "        df_plot.drop(['CD20','CD8','CD4','CK19','pH3'],axis=1,inplace=True)\n",
    "        #barplot\n",
    "        s_title = f'Relative {s_name}: Superpixel {k},{tu_level[0].replace(\"minimum\",\"\")}'.replace('Superpixel 1,','Std./CyCIF')\n",
    "        s_ylabel = f\"Std./CycIF Intensity\"\n",
    "        s_figname =f'{rootdir}/Figures/{s_measure}_Ratio_Single-Cyclic_0-3_manual_suppix{k}_{tu_level[0]}.png'\n",
    "        fig, ax = plt.subplots(figsize=(7,4),dpi=300)\n",
    "        sns.barplot(data=df_plot,ax=ax,palette='muted',order=ls_order_sp,ci='sd')\n",
    "        sns.stripplot(data=df_plot,ax=ax,palette='dark',order=ls_order_sp)\n",
    "        ax.set_ylabel(s_ylabel)\n",
    "        ax.set_xlabel(\"Marker\")\n",
    "        ax.set_ylim(0,3.5)\n",
    "        ax.axhline(1,color='black')\n",
    "        labels=ax.get_xticklabels()\n",
    "        ax.set_xticklabels(labels, rotation=90)\n",
    "        ax.set_title(s_title)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(s_figname)\n",
    "        if s_measure == 'mean_intensity_SBR':\n",
    "            ax.set_ylim(0,2.4)\n",
    "            ax.set_title(s_title) #,fontweight='bold',pad=10\n",
    "            fig.savefig(s_figname.replace('.png','_zoom.png'))\n",
    "        d_result.update({s_measure:df_result})\n",
    "        d_plot.update({s_measure:df_plot})\n",
    "        #break\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#alternative region\n",
    "ls_col = ['CD44', 'CD45', 'CK14', 'CK5', 'CK7', 'ER', 'Ecad','LaminAC', 'aSMA','Vim',] #good markers in norm breast\n",
    "for s_measure in ['mean_intensity','quartiles-0','quartiles-1','quartiles-2','quartiles-3','quartiles-4']:\n",
    "    df_alt.loc[:,f'{s_measure}_SBR'] = df_alt.loc[:,f'{s_measure}_fg']/df_alt.loc[:,'mean_intensity_bg']\n",
    "df_alt['dynamic_range_SBR'] = df_alt.loc[:,'quartiles-4_SBR'] - df_alt.loc[:,'quartiles-0_SBR']\n",
    "df_alt['dynamic_range'] = df_alt.loc[:,'quartiles-4_fg'] - df_alt.loc[:,'quartiles-0_fg'] #\n",
    "df_alt['tissue_label'] = df_alt.experimenttype + '_' + df_alt.label.astype('str')\n",
    "\n",
    "d_marker = {'normal':['CD44', 'CD45', 'CD68', 'CK14', 'CK5', 'CK7', 'ER', 'Ecad', 'Ki67',\n",
    "       'LaminAC', 'PCNA',  'Vim', 'aSMA'], #'PD1',\n",
    "            'tumor':['CD44', 'CD45', 'CD68',  'CK7',  'Ecad',  'Ki67', #'Her2',\n",
    "       'LaminAC', 'PCNA', 'PD1', 'Vim','aSMA' ]#'ER', #pd1 has some lamAC bleed through\n",
    "           }\n",
    "for s_type in ['normal','tumor']:\n",
    "    df_type = df_alt[df_alt.experimenttype == s_type]\n",
    "    ls_col = d_marker[s_type]\n",
    "    for s_measure in ['mean_intensity_fg','mean_intensity_bg']: #'dynamic_range','mean_intensity_SBR',\n",
    "        print(s_measure)\n",
    "        df_plot_alt = pd.DataFrame()\n",
    "        for s_marker in ls_col: #ls_marker:\n",
    "            df_marker = df_type[df_type.marker==s_marker]\n",
    "            df = df_marker[df_marker.experiment=='cyclic'].merge(df_marker[df_marker.experiment=='single'],on='tissue_label',suffixes=('_cyclic','_single'))\n",
    "            d_rename = dict(zip([i for i in range(len(df))],[item.split('_')[0] for item in df.tissue_label]))\n",
    "            se_sing = df.loc[:,f'{s_measure}_single'].rename(d_rename)\n",
    "            se_cyc = df.loc[:,f'{s_measure}_cyclic'].rename(d_rename)\n",
    "            df_plot_alt[s_marker] = (se_sing/se_cyc)\n",
    "        df_plot = d_plot[s_measure]        \n",
    "        df=pd.DataFrame()\n",
    "        df['ROI1'] = df_plot.loc[s_type,ls_col]\n",
    "        df['ROI2'] = df_plot_alt.loc[s_type,ls_col]\n",
    "        fig,ax=plt.subplots(figsize=(6,4),dpi=300)#\n",
    "        sns.scatterplot(x='ROI1', y='ROI2', data=df,hue=df.index,palette='tab20', ax=ax)\n",
    "        ax.legend(ncol=2,bbox_to_anchor=(1,1.05),fontsize=12,title='Marker')\n",
    "        #ax.set_ylim(0,2.4)\n",
    "        #ax.set_xlim(0,2.4)\n",
    "        r, pvalue = pearsonr(x=df.ROI1, y=df.ROI2)\n",
    "        ax.set_title(f'{d_measure[s_measure]} {s_type}\\n R = {r:.02}\\n p = {pvalue:.02}')\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'./Figures/Scatterplot_single_vs_cyclic_alt_{s_measure}_{s_type}.png')\n",
    "        print(r)\n",
    "        print(pvalue)\n",
    "        #break\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for s_measure, df_result in d_result.items(): \n",
    "    df_result.rename({0:'tumor',1:'normal'},axis=1,inplace=True)\n",
    "    try:\n",
    "        df_result.drop(['pH3_single','pH3_cyclic'],inplace=True)\n",
    "    except:\n",
    "        print('')\n",
    "    df = df_result.unstack().reset_index().dropna().rename({'level_0':'tissue','level_1':'marker_type',0:'SBR'},axis=1)\n",
    "    df['type'] = [item.split('_')[1] for item in df.marker_type]\n",
    "    df['marker'] = [item.split('_')[0] for item in df.marker_type]\n",
    "    df['marker_tissue'] = df.marker + '_' + df.tissue\n",
    "    #plot\n",
    "    df_plot = df[df.type=='single'].loc[:,['marker_tissue','SBR']].merge(df[df.type=='cyclic'].loc[:,['marker_tissue','SBR']],on='marker_tissue',suffixes=('_single','_cyclic'))\n",
    "    df_plot['marker_tissue'] = df_plot.marker_tissue.replace({'CK5_tumor':'CK5_normal'})\n",
    "    df_plot['marker'] = [item.split('_')[0] for item in df_plot.marker_tissue]\n",
    "    df_plot['tissue'] = [item.split('_')[1] for item in df_plot.marker_tissue]\n",
    "    df_plot.sort_values('marker_tissue',inplace=True)\n",
    "    df_plot.marker = df_plot.marker.replace({'Her2':'HER2'})\n",
    "    fig,ax=plt.subplots(figsize=(5.5,4),dpi=300)#\n",
    "    sns.scatterplot(x='SBR_cyclic', y='SBR_single', data=df_plot,hue='marker',palette='tab20', ax=ax)#,s=15\n",
    "    ax.legend(ncol=2,bbox_to_anchor=(1,1.05),fontsize=12,title='Marker')\n",
    "    r, pvalue = pearsonr(x=df_plot.SBR_cyclic, y=df_plot.SBR_single)\n",
    "    ax.set_ylabel('Standard IF')\n",
    "    ax.set_xlabel('CyCIF')\n",
    "    ax.set_ylim(0,22)\n",
    "    ax.set_xlim(0,22)\n",
    "    ax.set_title(f'{d_measure[s_measure]}\\n R = {r:.02}\\n p = {pvalue:.02}')\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Scatterplot_single_vs_cyclic_{s_measure}.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Percent Positive\n",
    "\n",
    "Thresholds were set on the images, using the pattern of staining to set an intensity value above which pixels were considered positive for a stain.\n",
    "\n",
    "Here, the same theresholds are applied to the single cell mean intensity values. If a cell has a mean intensity above the pixel threshold, it is considered positive. Percent of positive cells over all segmented cells is calculated for each marker."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#cropping used for images (different from the features)\n",
    "#new 20191209\n",
    "d_crop_image = {'44290-112':(5000,8800,1300,1800),\n",
    " '44290-113':(5886,8111,1300,1800),\n",
    " '44290-114':(5010,9242,1300,1800),\n",
    " '44290-115':(5975,10490,1300,1800),\n",
    " '44290-116':(6336,9174,1300,1800),\n",
    " '44294-116':(9547,5459,1300,1800),\n",
    " '44294-117':(12180,6092, 1300,1800),\n",
    " '44294-118':(10901, 6582, 1300,1800),\n",
    " '44294-119':(9853,5557, 1300,1800),\n",
    " '44294-120':(10250,5510, 1300,1800),\n",
    " }\n",
    "d_crop_feat = {'44290-112':(5000,8800,1300,1800),\n",
    " '44290-113':(2650,6400,1300,1800),\n",
    " '44290-114':(2650,6250,1300,1800),\n",
    " '44290-115':(2850,6200,1300,1800),\n",
    " '44290-116':(2650,6250,1300,1800),\n",
    " '44294-116':(7200,3580,12725,12365),\n",
    " '44294-117':(0,0,0,0),\n",
    " '44294-118':(0,0,0,0),\n",
    " '44294-119':(0,0,0,0),\n",
    " '44294-120':(0,0,0,0),\n",
    " }\n",
    "#define rois (see notes above about ROIs to exclude autofluorescence)\n",
    "df_rois = pd.DataFrame()\n",
    "for s_tissue, t_coord in d_crop_feat.items():\n",
    "    #print(s_tissue)\n",
    "    if t_coord == (0,0,0,0):\n",
    "        df_scene = df_mi[df_mi.slide==s_tissue]\n",
    "    else:\n",
    "        x_min = t_coord[0]\n",
    "        x_max = t_coord[0] + t_coord[2]\n",
    "        y_min = t_coord[1]\n",
    "        y_max = t_coord[1] + t_coord[3]\n",
    "        df_scene = df_mi[df_mi.slide==s_tissue]\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_scene)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate positive cells based on thresholding the mean intensity dataframe\n",
    "\n",
    "#empty dataframe\n",
    "df_pos = pd.DataFrame()\n",
    "\n",
    "#for each sample\n",
    "for s_index in df_t.index:\n",
    "    s_exp = df_t.loc[s_index,'experiment']\n",
    "    s_tissue = df_t.loc[s_index,'scene']\n",
    "    s_marker = df_t.loc[s_index,'marker']\n",
    "    s_type = df_t.loc[s_index,'experimenttype']\n",
    "    df_m_slide = df_m[(df_m.slide==s_tissue) & (df_m.marker == s_marker)]\n",
    "    if len(df_m_slide) !=1:\n",
    "        print(df_m_slide)\n",
    "    #when the features were extracted we had some exposure normalization ... need to undo\n",
    "    i_thresh = df_t.loc[s_index,'minimum']/(df_m_slide.loc[:,'refexp']/df_m_slide.loc[:,'exposure'])\n",
    "    s_marker_loc = df_rois.columns[pd.Series([item.split('_')[0]==s_marker for item in df_rois.columns])][0]\n",
    "    df_pos.loc[s_marker,f'{s_exp}_{s_type}'] = (df_rois[df_rois.slide==s_tissue].loc[:,s_marker_loc] > i_thresh[0]).sum()\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#generate dataframe for plotting\n",
    "#tumor\n",
    "df_percent = pd.DataFrame()\n",
    "df_percent['tumor'] = (df_pos.dropna(axis=0).single_tumor/df_pos.dropna(axis=0).cyclic_tumor).fillna(1)\n",
    "#normal\n",
    "df_percent['normal'] = (df_pos.dropna(axis=0).single_normal/df_pos.dropna(axis=0).cyclic_normal).fillna(1)\n",
    "\n",
    "\n",
    "\n",
    "df_tum = pd.DataFrame(df_percent.unstack().tumor, columns=['value'])\n",
    "df_norm = pd.DataFrame(df_percent.unstack().normal, columns=['value'])\n",
    "df_tum['type'] = 'tumor'\n",
    "df_norm['type'] = 'normal'\n",
    "df_tum['marker'] = df_tum.index.tolist()\n",
    "df_norm['marker'] = df_norm.index.tolist()\n",
    "df_longer = df_tum.append(df_norm)\n",
    "df_longer = df_longer.drop('pH3')\n",
    "\n",
    "ls_order = df_tum.index\n",
    "\n",
    "#stats\n",
    "print(df_longer.groupby('marker').mean().mean())\n",
    "print(df_longer.groupby('marker').mean().std())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#barplot\n",
    "\n",
    "s_title = 'Relative % Pos. of Standard versus CycIF'\n",
    "s_ylabel = \"Std./CycIF % Pos.\"\n",
    "s_figname =f'{codedir}/Figures/Percent_Positive_Ratio_Single-Cyclic.png'\n",
    "fig, ax = plt.subplots(figsize=(7,4),dpi=300)\n",
    "sns.barplot(data=df_longer, x='marker', y='value', palette = 'muted',ci='sd')\n",
    "sns.stripplot(data=df_longer, x='marker', y='value', palette = 'dark')\n",
    "ax.set_ylabel(s_ylabel)\n",
    "\n",
    "ax.set_xlabel(\"Marker\")\n",
    "ax.axhline(1,color='black')\n",
    "\n",
    "labels=ax.get_xticklabels()\n",
    "ax.set_xticklabels(labels, rotation=90)\n",
    "ax.set_title(s_title)\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.savefig(s_figname)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df = df_result.unstack().reset_index().dropna().rename({'level_0':'tissue','level_1':'marker_type',0:'SBR'},axis=1)\n",
    "df = df_pos.unstack().reset_index().dropna().rename({'level_0':'type_tissue','level_1':'marker',0:'Positive'},axis=1)\n",
    "df = df[df.marker!='pH3']   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['type'] = [item.split('_')[0] for item in df.type_tissue]\n",
    "df['tissue'] = [item.split('_')[1] for item in df.type_tissue]\n",
    "df['marker_tissue'] = df.marker + '_' + df.tissue\n",
    "df_plot = df[df.type=='single'].loc[:,['marker_tissue','Positive']].merge(df[df.type=='cyclic'].loc[:,['marker_tissue','Positive']],on='marker_tissue',suffixes=('_single','_cyclic'))\n",
    "df_plot['marker_tissue'] = df_plot.marker_tissue.replace({'CK5_tumor':'CK5_normal'})  \n",
    "df_plot['marker'] = [item.split('_')[0] for item in df_plot.marker_tissue]\n",
    "df_plot['tissue'] = [item.split('_')[1] for item in df_plot.marker_tissue]\n",
    "df_plot.sort_values('marker_tissue',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(3,4),dpi=300)\n",
    "sns.scatterplot(x='Positive_cyclic', y='Positive_single', data=df_plot,hue='marker',palette='muted', ax=ax)\n",
    "ax.legend(ncol=2,bbox_to_anchor=(1,1.05),fontsize=10,title='Marker')\n",
    "ax.get_legend().remove()\n",
    "r, pvalue = pearsonr(x=df_plot.Positive_cyclic, y=df_plot.Positive_single)\n",
    "ax.set_ylabel('Standard IF')\n",
    "ax.set_xlabel('CyCIF')\n",
    "ax.ticklabel_format(axis='both',style='sci',scilimits=(0,0))\n",
    "ax.set_title(f'Number Positive \\n R = {r:.02}\\n p = {pvalue:.02}',pad=20)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./Figures/Scatterplot_single_vs_cyclic_Positive.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(df.Positive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate percent positive cells based on thresholding the mean intensity dataframe\n",
    "\n",
    "#empty dataframe\n",
    "df_pos = pd.DataFrame()\n",
    "\n",
    "#for each sample\n",
    "for s_index in df_t.index:\n",
    "    s_exp = df_t.loc[s_index,'experiment']\n",
    "    s_tissue = df_t.loc[s_index,'scene']\n",
    "    s_marker = df_t.loc[s_index,'marker']\n",
    "    s_type = df_t.loc[s_index,'experimenttype']\n",
    "    df_m_slide = df_m[(df_m.slide==s_tissue) & (df_m.marker == s_marker)]\n",
    "    if len(df_m_slide) !=1:\n",
    "        print(df_m_slide)\n",
    "    #when the features were extracted we had some exposure normalization ... need to undo\n",
    "    i_thresh = df_t.loc[s_index,'minimum']/(df_m_slide.loc[:,'refexp']/df_m_slide.loc[:,'exposure'])\n",
    "    s_marker_loc = df_rois.columns[pd.Series([item.split('_')[0]==s_marker for item in df_rois.columns])][0]\n",
    "    df_pos.loc[s_marker,f'{s_exp}_{s_type}'] = (df_rois[df_rois.slide==s_tissue].loc[:,s_marker_loc] > i_thresh[0]).sum()/len(df_rois[df_rois.slide==s_tissue])\n",
    "    #break\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df = df_result.unstack().reset_index().dropna().rename({'level_0':'tissue','level_1':'marker_type',0:'SBR'},axis=1)\n",
    "df = df_pos.unstack().reset_index().dropna().rename({'level_0':'type_tissue','level_1':'marker',0:'Positive'},axis=1)\n",
    "df = df[df.marker!='pH3']   \n",
    "df['type'] = [item.split('_')[0] for item in df.type_tissue]\n",
    "df['tissue'] = [item.split('_')[1] for item in df.type_tissue]\n",
    "df['marker_tissue'] = df.marker + '_' + df.tissue\n",
    "df_plot = df[df.type=='single'].loc[:,['marker_tissue','Positive']].merge(df[df.type=='cyclic'].loc[:,['marker_tissue','Positive']],on='marker_tissue',suffixes=('_single','_cyclic'))\n",
    "df_plot['marker_tissue'] = df_plot.marker_tissue.replace({'CK5_tumor':'CK5_normal'})  \n",
    "df_plot['marker'] = [item.split('_')[0] for item in df_plot.marker_tissue]\n",
    "df_plot['tissue'] = [item.split('_')[1] for item in df_plot.marker_tissue]\n",
    "df_plot.sort_values('marker_tissue',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(3.2,4),dpi=300)\n",
    "sns.scatterplot(x='Positive_cyclic', y='Positive_single', data=df_plot,hue='marker',palette='muted', ax=ax)\n",
    "ax.legend(ncol=2,bbox_to_anchor=(1,1.05),fontsize=10,title='Marker')\n",
    "ax.get_legend().remove()\n",
    "r, pvalue = pearsonr(x=df_plot.Positive_cyclic, y=df_plot.Positive_single)\n",
    "ax.set_ylabel('Standard IF')\n",
    "ax.set_xlabel('CyCIF')\n",
    "#ax.ticklabel_format(axis='both',style='sci',scilimits=(0,0))\n",
    "ax.set_title(f'Fraction Positive \\n R = {r:.02}\\n p = {pvalue:.02}',pad=20)\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_ylim(0,1)\n",
    "plt.tight_layout()\n",
    "#fig.savefig(f'./Figures/Scatterplot_single_vs_cyclic_Positive.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.histplot(df.Positive.rename('Fraction Positive'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tissue loss <a name=\"tissue\"></a>\n",
    "\n",
    "[contents](#contents)\n",
    "\n",
    "quantify tissue loss across 10 rounds of quenching, normal/tumor tissue: Biomax FDA808l-2 : Multiple tumor (24 organs )and normal (6 organs) tissue array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load tissue loss data\n",
    "\n",
    "s_slides = '808L2'# 'HER2A'#\n",
    "s_in = f'{codedir}/Data/Quench/features_{s_slides}_TissueLoss.csv'\n",
    "df = pd.read_csv(s_in,index_col=0,low_memory=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#add scene\n",
    "df['scene'] = [item.split('_')[1] for item in df.slide_scene_new]\n",
    "#tissue completely lost\n",
    "ls_na = ['DAPI10_nuclei_thresh','DAPI7_nuclei_thresh','DAPI8_nuclei_thresh','DAPI9_nuclei_thresh']\n",
    "df.loc[:,ls_na] = df.loc[:,ls_na].fillna(False)\n",
    "for s_col in ['nuclei_area', 'nuclei_eccentricity']:\n",
    "    df[f'{s_col}_Q'] = pd.qcut(df.loc[:,s_col],4,labels=[\"Q1\", \"Q2\", \"Q3\",\"Q4\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load annotation\n",
    "df_a = pd.read_csv(f'{codedir}/Data/FDA808l-2 specs.csv')\n",
    "#grade\n",
    "df_a['Grade'] = df_a.Grade.replace('-',np.nan)\n",
    "#stage\n",
    "d_replace = {'-':np.nan, 'IA':'I', 'III':\"III\", 'IIA':\"II\", 'IIB':'II', 'IIIA':\"III\", 'IB':\"I\",\n",
    "       ' IIB G2':\"II\", 'IIIB':\"III\", 'IVB':'IV', 'IIA G3':\"II\", 'IIB G3 ':\"II\", 'IIB G2':\"II\"}\n",
    "df_a['Stage'] = df_a.Stage.replace(d_replace)\n",
    "df_a.loc[:,'Pathology'] = df_a.loc[:,'Pathology diagnosis']\n",
    "df_a.loc[:,'Age_Q'] = pd.qcut(df_a.Age,4,labels=[\"Q1\", \"Q2\", \"Q3\",\"Q4\"])\n",
    "df_a.loc[~(df_a.Source.str.contains('\\*')),'Source'] = 'autopsy'\n",
    "df_a.loc[df_a.Source.str.contains('\\*'),'Source'] = 'resection'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#shorten names\n",
    "d_rename = dict(zip(df.columns[df.dtypes=='bool'],[item.split('_')[0] for item in df.columns[df.dtypes=='bool']]))\n",
    "df.rename(d_rename,axis=1,inplace=True)\n",
    "ls_order = [ 'DAPI1', 'DAPI2', 'DAPI3', 'DAPI4',\n",
    "       'DAPI5', 'DAPI6', 'DAPI7', 'DAPI8', 'DAPI9','DAPI10'] #'DAPI-1', 'DAPI0',"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate fraction per scene, round\n",
    "ls_drop = ['scene006','scene012','scene013','scene016', 'scene072'] #scene6 lost in 55, 72 mis registered in 53, others few cells\n",
    "df_all = pd.DataFrame()\n",
    "fig, ax = plt.subplots(figsize=(8,4),dpi=300)\n",
    "for s_scene in sorted(set(df.scene.unique()) - set(ls_drop)):\n",
    "    df_sum = df[df.scene==s_scene].loc[:,ls_order+['slide']].groupby('slide').sum()\n",
    "    se_max = df_sum.max(axis=1) #df_sum.loc[:,'DAPI1'] #\n",
    "    df_plot = (df_sum.T/se_max) #normalize\n",
    "    #QC\n",
    "    if len (df_plot.std().loc[df_plot.std() > .1]) > 1:\n",
    "        print(s_scene)\n",
    "        print(df_plot.std().loc[df_plot.std() > .2].index)\n",
    "    if sum(df_plot.loc['DAPI3'] - df_plot.loc['DAPI2'] > .1) > 0:\n",
    "        print(s_scene)\n",
    "    df_long = pd.DataFrame(df_plot.unstack()).reset_index().rename({'level_1':'Round',0:'Fraction_Cells'},axis=1)\n",
    "    sns.lineplot(x='Round',y='Fraction_Cells',data=df_long,ax=ax)\n",
    "    df_long['scene'] = s_scene\n",
    "    df_all = df_all.append(df_long)\n",
    "    #break\n",
    "ax.set_ylim(0,1)\n",
    "ax.set_xticks(ls_order)\n",
    "ax.set_xticklabels(ls_order, rotation=90)\n",
    "ax.set_title('Tissue Retention per Round, per Core')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./Figures/Tissue_Loss_per_Scene_all.png')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#all \n",
    "fig, ax = plt.subplots(figsize=(5,3),dpi=300)\n",
    "#sns.stripplot(x='Round',y='Fraction_Cells',data=df_all.reset_index(),ax=ax,color='C0',s=1,alpha=0.3,jitter=.2)\n",
    "sns.lineplot(x='Round',y='Fraction_Cells',data=df_all.reset_index(),ax=ax) #,estimator='scene'\n",
    "ax.set_ylim(.8,1)\n",
    "ax.set_xticks(ls_order)\n",
    "ax.set_xticklabels([item.replace('DAPI','') for item in ls_order])\n",
    "ax.set_ylabel('Fraction Cells')\n",
    "ax.set_title('Tissue Retention by Round') #,fontweight='bold'\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./Figures/Tissue_Loss_per_Round.png')    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#cats\n",
    "d_result = {}\n",
    "for s_type in ['Type','Grade','Stage','Sex','Organ','Pathology','Age_Q']: #\n",
    "    fig, ax = plt.subplots(figsize=(6,3),dpi=300)\n",
    "    df_all['Cat'] = df_all.scene.map(dict(zip(df_a.Scene,df_a.loc[:,s_type])))\n",
    "    sns.lineplot(x='Round',y='Fraction_Cells',hue='Cat',data=df_all.reset_index(),ax=ax) #,estimator='scene'\n",
    "    ax.set_xticks(ls_order)\n",
    "    ax.set_xticklabels([item.replace('DAPI','') for item in ls_order])\n",
    "    ax.set_ylabel('Fraction Cells')\n",
    "    ax.set_ylim(.75,1)\n",
    "    ax.set_title(f'Tissue Retention by {s_type}',pad=10) #fontweight='bold',\n",
    "    ax.legend(bbox_to_anchor=(.95,.95))\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Tissue_Loss_per_Round_{s_type}.png')    \n",
    "    #test\n",
    "    df_last = df_all[df_all.Round=='DAPI10'].copy()\n",
    "    df_last.index = df_all[df_all.Round=='DAPI10'].slide + '_' + df_all[df_all.Round=='DAPI10'].scene\n",
    "    df_first = df_all[df_all.Round=='DAPI1'].copy()\n",
    "    df_first.index = df_all[df_all.Round=='DAPI1'].slide + '_' + df_all[df_all.Round=='DAPI1'].scene\n",
    "    df_first['Remaining'] = df_last.loc[:,'Fraction_Cells']/df_first.loc[:,'Fraction_Cells']\n",
    "    figsize = (3.5,3)\n",
    "    if s_type == 'Type':\n",
    "        figsize = (4.8,3)\n",
    "    fig,ax=plt.subplots(figsize=figsize,dpi=300)\n",
    "    sns.boxplot(x='Cat',y='Remaining',data=df_first,order = df_first.Cat.dropna().sort_values().unique(),fliersize=0,ax=ax)\n",
    "    sns.stripplot(x='Cat',y='Remaining',data=df_first,order = df_first.Cat.dropna().sort_values().unique(),palette='dark',ax=ax)\n",
    "    #anova\n",
    "    lls_f = [df_first.Remaining[df_first.Cat == item] for item in df_first.Cat.dropna().unique()]\n",
    "    #statistic,pvalue = stats.f_oneway(*lls_f)\n",
    "    statistic,pvalue = stats.kruskal(*lls_f)\n",
    "    ax.set_ylabel('Fraction Cells \\n After 10 Rounds')\n",
    "    ax.set_title(f'{s_type} \\n p = {pvalue:.02}') #pad=10,,fontweight='bold' Ten Round Retention by \n",
    "    ax.set_xlabel(s_type)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Tissue_Loss_10round_{s_type}.png') \n",
    "    d_result.update({s_type:df_first})\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#surgical versus autopsy\n",
    "s_type = 'Source'\n",
    "df_first = d_result['Type'][d_result['Type'].Cat=='Normal'].copy()\n",
    "df_first['Cat'] = df_first.scene.map(dict(zip(df_a.Scene,df_a.Source)))\n",
    "#test\n",
    "df_first['Remaining'] = df_last.loc[:,'Fraction_Cells']/df_first.loc[:,'Fraction_Cells']\n",
    "fig,ax=plt.subplots(figsize=(3.6,3),dpi=300)\n",
    "sns.boxplot(x='Cat',y='Remaining',data=df_first,order = df_first.Cat.dropna().sort_values().unique(),fliersize=0,ax=ax)\n",
    "sns.stripplot(x='Cat',y='Remaining',data=df_first,order = df_first.Cat.dropna().sort_values().unique(),palette='dark',ax=ax)\n",
    "#anova\n",
    "lls_f = [df_first.Remaining[df_first.Cat == item] for item in df_first.Cat.dropna().unique()]\n",
    "statistic,pvalue = stats.mannwhitneyu(*lls_f)\n",
    "#statistic,pvalue = stats.f_oneway(*lls_f)\n",
    "ax.set_ylabel('Fraction Cells \\n After 10 Rounds')\n",
    "ax.set_title(f'{s_type} \\n p = {pvalue:.02}') #pad=10,,fontweight='bold' Ten Round Retention by \n",
    "ax.set_xlabel(s_type)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./Figures/Tissue_Loss_10round_{s_type}.png') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#treands between normal and malignant\n",
    "for s_type in ['Sex','Organ','Pathology','Age_Q','Grade','Type','Stage']:#,\n",
    "    df_first = d_result[s_type]\n",
    "    df_first['Type'] = df_first.scene.map(dict(zip(df_a.Scene,df_a.Type)))\n",
    "    for s_tum in ['Normal','Malignant']:\n",
    "        df_plot = df_first[df_first.Type==s_tum]\n",
    "        if len(df_plot.dropna().Cat.unique()) > 1:\n",
    "            width = 6\n",
    "            if s_tum=='Malignant':\n",
    "                if s_type == 'Pathology':\n",
    "                    width = 8\n",
    "            fig,ax=plt.subplots(figsize=(width,len(df_plot.Cat.unique())*.3+1.5))\n",
    "            order =df_plot.groupby('Cat').mean().sort_values(by='Remaining').index\n",
    "            sns.stripplot(y='Cat',x='Remaining',data=df_plot,order=order,palette='muted',ax=ax)\n",
    "            sns.pointplot(y='Cat',x='Remaining',data=df_plot,order=order,ax=ax,\n",
    "                  palette=\"dark\",markers=\"d\", scale=1, ci=95)\n",
    "            #anova\n",
    "            lls_f = [df_plot.Remaining[df_plot.Cat == item] for item in df_plot.Cat.dropna().unique()]\n",
    "            #statistic,pvalue = stats.f_oneway(*lls_f)\n",
    "            statistic,pvalue = stats.kruskal(*lls_f)\n",
    "            ax.set_title(f'{s_tum} Tissue \\n p = {pvalue:.02}')\n",
    "            ax.set_ylabel(s_type)\n",
    "            ax.set_xlim(0,1.1)\n",
    "            plt.tight_layout()\n",
    "            fig.savefig(f'./Figures/Tissue_Loss_10round_{s_type}_{s_tum}.png') \n",
    "            #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate fraction per scene, round, cell size or shape\n",
    "for s_col in ['nuclei_area', 'nuclei_eccentricity']:\n",
    "    df[f'{s_col}_Q'] = pd.qcut(df.loc[:,s_col],4,labels=[\"Q1\", \"Q2\", \"Q3\",\"Q4\"])\n",
    "ls_drop = ['scene006','scene012','scene013','scene016', 'scene072'] #scene6 lost in 55, 72 mis registered in 53, others few cells\n",
    "df_cell = pd.DataFrame()\n",
    "\n",
    "for s_scene in sorted(set(df.scene.unique()) - set(ls_drop)):\n",
    "    #print(s_scene)\n",
    "    for s_cell in ['nuclei_area_Q','nuclei_eccentricity_Q']:\n",
    "        #print(s_cell)\n",
    "        df_sum = df[df.scene==s_scene].loc[:,ls_order+['slide',s_cell]].groupby(['slide',s_cell]).sum()\n",
    "        df_max = pd.DataFrame(df_sum.reset_index().groupby(['slide',s_cell]).max().max(axis=1))\n",
    "        for s_order in ls_order:\n",
    "            df_max.loc[:,s_order] = df_max.loc[:,0]\n",
    "        df_max.drop(0,axis=1,inplace=True)\n",
    "        df_plot = (df_sum/df_max) #normalize\n",
    "        df_long = df_plot.unstack().unstack().reset_index().rename({'level_0':'Round',0:'Fraction_Cells',s_cell:'Quartile'},axis=1)\n",
    "        df_long['scene'] = s_scene\n",
    "        df_long['cell'] = s_cell\n",
    "        df_cell = df_cell.append(df_long)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cats\n",
    "d_result_cell = {}\n",
    "for s_type in ['nuclei_area_Q','nuclei_eccentricity_Q']: #\n",
    "    fig, ax = plt.subplots(figsize=(4.5,3),dpi=300)\n",
    "    #df_cell['Cat'] = df_all.scene.map(dict(zip(df_a.Scene,df_a.loc[:,s_type])))\n",
    "    sns.lineplot(x='Round',y='Fraction_Cells',hue='Quartile',data=df_cell[df_cell.cell==s_type].reset_index(),ax=ax) #\n",
    "    ax.set_xticks(ls_order)\n",
    "    ax.set_xticklabels([item.replace('DAPI','') for item in ls_order])\n",
    "    ax.set_ylabel('Fraction Cells')\n",
    "    ax.set_ylim(.8,1)\n",
    "    ax.set_title(f'Tissue Retention by {s_type.replace(\"_\",\" \").replace(\"Q\",\"\")}',pad=10) #fontweight='bold',\n",
    "    ax.legend(bbox_to_anchor=(.95,.95))\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Tissue_Loss_per_Round_{s_type}.png')    \n",
    "\n",
    "    #test\n",
    "    df_cel = df_cell[df_cell.cell==s_type]\n",
    "    df_last = df_cel[df_cel.Round=='DAPI10'].copy()\n",
    "    df_last.index = df_cel[df_cel.Round=='DAPI10'].slide + '_' + df_cel[df_cel.Round=='DAPI10'].scene\n",
    "    df_first = df_cel[df_cel.Round=='DAPI1'].copy()\n",
    "    df_first.index = df_cel[df_cel.Round=='DAPI1'].slide + '_' + df_cel[df_cel.Round=='DAPI1'].scene\n",
    "    df_first['Remaining'] = df_last.loc[:,'Fraction_Cells']/df_first.loc[:,'Fraction_Cells']\n",
    "    fig,ax=plt.subplots(figsize=(4,3),dpi=200)\n",
    "    sns.boxplot(x='Quartile',y='Remaining',data=df_first,order = df_first.Quartile.dropna().sort_values().unique(),fliersize=0,ax=ax)\n",
    "    sns.stripplot(x='Quartile',y='Remaining',data=df_first,order = df_first.Quartile.dropna().sort_values().unique(),palette='dark',ax=ax)\n",
    "    #anova\n",
    "    lls_f = [df_first.Remaining[df_first.Quartile == item] for item in df_first.Quartile.dropna().unique()]\n",
    "    #statistic,pvalue = stats.f_oneway(*lls_f)\n",
    "    statistic,pvalue = stats.kruskal(*lls_f)\n",
    "    ax.set_ylabel('Fraction Cells \\n After 10 Rounds')\n",
    "    #ax.set_title(f'Ten Round Retention by {s_type} \\n p = {pvalue:.02}',fontweight='bold') #pad=10,\n",
    "    ax.set_title(f'{s_type.replace(\"_\",\" \").replace(\"Q\",\"\").replace(\"nuc\",\" Nuc\")}\\n p = {pvalue:.02}')#Ten Round Retention \\n by,pad=10 ,fontweight='bold'\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Tissue_Loss_10round_{s_type}.png') \n",
    "    d_result_cell.update({s_type:df_first})\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#treands between normal and malignant\n",
    "for s_type in ['nuclei_area_Q','nuclei_eccentricity_Q']:#,\n",
    "    df_first = d_result_cell[s_type]\n",
    "    df_first['Type'] = df_first.scene.map(dict(zip(df_a.Scene,df_a.Type)))\n",
    "    for s_tum in ['Normal','Malignant']:\n",
    "        df_plot = df_first[df_first.Type==s_tum]\n",
    "        if len(df_plot.dropna().Quartile.unique()) > 1:\n",
    "            fig,ax=plt.subplots(figsize=(5,len(df_plot.Quartile.unique())*.3+1.5),dpi=300)\n",
    "            order =df_plot.groupby('Quartile').mean().sort_values(by='Remaining').index\n",
    "            sns.stripplot(y='Quartile',x='Remaining',data=df_plot,order=order,palette='muted',ax=ax)\n",
    "            sns.pointplot(y='Quartile',x='Remaining',data=df_plot,order=order,ax=ax,\n",
    "                  palette=\"dark\",markers=\"d\", scale=1, ci=95)\n",
    "            #anova\n",
    "            lls_f = [df_plot.Remaining[df_plot.Quartile == item] for item in df_plot.Quartile.dropna().unique()]\n",
    "            #statistic,pvalue = stats.f_oneway(*lls_f)\n",
    "            statistic,pvalue = stats.kruskal(*lls_f)\n",
    "            ax.set_title(f'{s_tum} Tissue \\n p = {pvalue:.02}')\n",
    "            ax.set_ylabel(s_type)\n",
    "            ax.set_xlim(0,1.1)\n",
    "            plt.tight_layout()\n",
    "            fig.savefig(f'./Figures/Tissue_Loss_10round_{s_type}_{s_tum}.png') \n",
    "            #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3.9.5",
   "language": "python",
   "name": "base"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}