Quenching_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of Quenching\n",
    "\n",
    "**Questions:** \n",
    "1. What is the optimal time, H2O2 concentration, and heat/light conditon for efficient signal removal.\n",
    "2. How does autofluorescence change over repeated quenching reactions?\n",
    "3. Is autofluorescence decrease due to repeated imaging, or quenching itself?\n",
    "\n",
    "**Samples:** \n",
    "\n",
    "***Experiment 1***\n",
    "\n",
    "Six adjacent sections of normal pancreas were processed with the cyclif IF tissue processing protocol. Three were stained with an antibody cocktail, and three were left as no-antbody controls. Slides we repeatedly quenched for 15 minutes and imaged. Quenching solution contained 3%, 4.5% and 6% H2O2.\n",
    "\n",
    "***Experiment 2***\n",
    "\n",
    "Four adjacent sections of breast cancer were processed, three were stained with strong antibody cocktail; one was left as no-antibody control. Stained slides were repeatedly quenched in 3% H202 for 30 minutes, one hour, or 30 minutes plus heat, respectively. Control was repeatedly quenched 30 minutes without heat. All slides imaged between every quenching interval. \n",
    "\n",
    "***Experiment 3***\n",
    "\n",
    "Two adjacent sections of normal pancreas were processed and left without antibody. They were imaged before any quenching was applied. Then one was repeatedly quenched in 3% H2O2 for 30 miutes with heat, while the other was only imaged.\n",
    "\n",
    "\n",
    "**Method**: Slides were imaged as described above, registered and tissue area was segmented out based on flooding DAPI mask. Mean intensity within the tissue area was quantififed.\n",
    "\n",
    "***Experiment 4***\n",
    "\n",
    "Single-cell clustering results from [this notebook](https://github.com/engjen/cycIF_Validation/blob/master/Fig2_Quenching_Single_Cell.ipynb) were imported for analysis of distinct cell populations' autofluorescence change over time.\n",
    "\n",
    "We analyzed 3 datasets:\n",
    "\n",
    "- Normal Pancreas Tissue (n=3, 4 rounds)\n",
    "\n",
    "- Multiple tumor (24 organs) and normal (6 organs) tissue array, with grade and TNM/Stage,72 cases/72 cores(core size 1.5mm) (n=216, 10 rounds)\n",
    "\n",
    "- Her-2 IHC Bio-marker control TMA, with breast invasive ductal carcinoma, including TNM, clinical stage and pathology grade, 10 cases/ 10 cores (n=30, 6 rounds)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib\n",
    "import seaborn as sns\n",
    "matplotlib.rcParams.update({'font.size': 15})\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "from scipy import stats"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "# Compare Single Cell Clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import single cell clustering\n",
    "#HER2A_AF_raw_20_leiden0.2.csv\n",
    "s_sample = 'HER2A_c2_nuclei'\n",
    "s_sample = '808L2_c2_nuclei'\n",
    "s_slide =  s_sample.split('_')[0] #'HER2A'#'808L2'#'4165NPanc'#\n",
    "s_slide = '4165NPanc'#\n",
    "n_neighbors= 25 #25 or 20\n",
    "resolution=0.2\n",
    "if s_slide == '4165NPanc':\n",
    "    df_cluster = pd.read_csv(f'{codedir}/Data/{s_slide}_AF_raw_leiden.csv',index_col=0)\n",
    "else:\n",
    "    df_cluster = pd.read_csv(f'{codedir}/Data/{s_sample}_raw_{n_neighbors}_leiden{resolution}.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#check\n",
    "df_cluster['scene'] = [item.split('_')[1] for item in df_cluster.index]\n",
    "df_cluster.sort_index(axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rename\n",
    "d_replace = {'R00c2_nuclei':0,'R0c2_nuclei':30,  'R1c2_nuclei':60, 'R2c2_nuclei':90, 'R3c2_nuclei':120,\n",
    "       'R4c2_nuclei':150, 'R5c2_nuclei':180, 'R6c2_nuclei':210,'R7c2_nuclei':240, 'R8c2_nuclei':270,\n",
    "             'R9c2_nuclei':300,'R10c2_nuclei':330,'R00c3_nuclei': 0, 'R0c3_nuclei': 30, 'R1c3_nuclei': 60, 'R2c3_nuclei': 90,\n",
    " 'R3c3_nuclei': 120, 'R4c3_nuclei': 150, 'R5c3_nuclei': 180, 'R6c3_nuclei': 210, 'R7c3_nuclei': 240, 'R8c3_nuclei': 270,\n",
    " 'R9c3_nuclei': 300, 'R10c3_nuclei': 330\n",
    "            }\n",
    "df_cluster.rename(d_replace,axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_order = df_cluster.columns[df_cluster.dtypes=='float64'].sort_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " df.groupby(['leiden']).mean().drop(ls_drop).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot clusters\n",
    "d_rename = {'c2':'AF488','c3':'A555'}\n",
    "matplotlib.rcParams.update({'font.size': 15})\n",
    "if s_slide == '4165NPanc':\n",
    "    d_color = {'nuc':'Nuclei','cyt':'Cytoplasm'}\n",
    "    ls_drop = [8,9,10]\n",
    "    ls_labels = ['0','15','30','45','60']\n",
    "    for s_loc in sorted(d_color.keys()):\n",
    "        s_nloc = sorted(set(d_color.keys()) - {s_loc})[0]\n",
    "        df = df_cluster.loc[:,~df_cluster.columns.str.contains(s_nloc)]\n",
    "        df_std = df.groupby(['leiden','tissue']).mean().drop(ls_drop)\n",
    "        df_std = df_std.groupby('leiden').std().T\n",
    "        fig,ax=plt.subplots(figsize=(6,4),dpi=300)\n",
    "        df.groupby(['leiden']).mean().drop(ls_drop).T.plot(yerr=df_std,ax=ax)\n",
    "        ax.set_xticks(ax.get_xticks()[1:-1])\n",
    "        ax.set_xticklabels(ls_labels,rotation=25)\n",
    "        ax.set_xlabel('Minutes of Quenching')\n",
    "        ax.set_ylabel('Mean Intensity')\n",
    "        ax.set_title(f'AF488 Autofluorescence \\nin {d_color[s_loc]}',fontsize=18)\n",
    "        ax.set_ylim(0,2200)\n",
    "        ax.legend(loc='upper left', bbox_to_anchor=(1.1, .95),ncol=1)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/Leiden_{s_slide}_{s_loc}_quenching.png')\n",
    "\n",
    "else:\n",
    "    if s_slide=='HER2A':\n",
    "        if s_sample == 'HER2A_c3_nuclei':\n",
    "            ls_drop = [8,9,10,11]\n",
    "        elif s_sample == 'HER2A_c2_nuclei':   \n",
    "            ls_drop = [8,9]\n",
    "        ylim = 0,2500\n",
    "        ncol=1\n",
    "        figsize = (5,4)\n",
    "        ls_scene = ['scene006','scene004','scene005','scene009']\n",
    "        df_cluster = df_cluster[df_cluster.scene.isin(ls_scene)]\n",
    "        #ls_labels = ['0','30','60','90','120','150'] #,'180'\n",
    "    if s_slide == '808L2':\n",
    "        if s_sample == '808L2_c2_nuclei':\n",
    "            ls_drop = [18,19]\n",
    "            ylim = 0,10000\n",
    "        elif s_sample == '808L2_c3_nuclei':\n",
    "            ls_drop = [15,16,17,18,20]\n",
    "        ncol = 2\n",
    "        figsize = (6,4)\n",
    "    ls_order = df_cluster.columns[df_cluster.dtypes=='float64'].sort_values()\n",
    "    df = df_cluster.drop('scene',axis=1)\n",
    "    df_std = df.groupby(['leiden','tissue']).mean().drop(ls_drop)\n",
    "    df_std = df_std.groupby('leiden').std().T\n",
    "    fig,ax=plt.subplots(figsize=figsize,dpi=300)\n",
    "    df.groupby(['leiden']).mean().drop(ls_drop).T.loc[ls_order].plot(yerr=df_std.loc[ls_order],ax=ax,colormap='tab20')\n",
    "    ax.set_xlabel('Minutes of Quenching')\n",
    "    ax.set_ylabel('Mean Intensity')\n",
    "    ax.set_title(f'{d_rename[s_sample.split(\"_\")[1]]} Autofluorescence \\nin {s_sample.split(\"_\")[2]}',fontsize=18)\n",
    "    ax.set_ylim(ylim)\n",
    "    ax.legend(loc='upper left', bbox_to_anchor=(.98, 1.01),ncol=ncol,fontsize=12) #\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{codedir}/Figures/Leiden_{s_sample}__quenching.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#120 minutes is lowest AF, 270 minutes is highest, 240 near 270, #180 is R5Q; in c3 240 is higher\n",
    "d_rename = {'c2':'AF488','c3':'A555'}\n",
    "if s_slide == '808L2':\n",
    "    i_low = 120\n",
    "    i_high = 240\n",
    "if s_slide == 'HER2A':\n",
    "    i_low = 90\n",
    "    i_high = 210\n",
    "matplotlib.rcParams.update({'font.size': 15})\n",
    "df_m =  df.groupby(['leiden','tissue']).mean().drop(ls_drop)\n",
    "df_diff = pd.DataFrame()\n",
    "ls_lei= sorted(set(df.leiden.unique()) - set(ls_drop))#[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]\n",
    "df.loc[:,df.dtypes=='float64'].mean().sort_values()\n",
    "for i_lei in ls_lei: \n",
    "    se_diff = (df_m.loc[i_lei,i_high] - df_m.loc[i_lei,i_low])/df_m.loc[i_lei,i_high]\n",
    "    se_diff.name= i_lei\n",
    "    df_diff = df_diff.append(se_diff)\n",
    "ls_order = df_diff.T.mean().sort_values(ascending=False).index.tolist()\n",
    "fig, ax = plt.subplots(figsize=(6,4),dpi=300)\n",
    "sns.barplot(data=df_diff.T,order=ls_order,ax=ax,palette='muted')\n",
    "sns.stripplot(data=df_diff.T,order=ls_order,ax=ax,palette='dark')\n",
    "ax.set_title(f'{d_rename[s_sample.split(\"_\")[1]]} Autofluorescence Remaining \\nAfter subtraction of {i_low} min AF')\n",
    "ax.set_ylabel(f'Fraction {i_high} min AF')\n",
    "ax.set_xlabel('Cluster')\n",
    "plt.tight_layout()\n",
    "print(1 - df_diff.mean().mean())\n",
    "fig.savefig(f'./Figures/Optimal_subtraction_{s_sample}.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## compare leiden to annotation\n",
    "\n",
    "808L2 has biomax annotation info"
   ]
  },
  {
   "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",
    "# add scene\n",
    "#adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#chi squared and post hoc\n",
    "matplotlib.rcParams.update({'font.size': 7})\n",
    "#ls drop from above\n",
    "ls_index = df_cluster[~df_cluster.leiden.isin(ls_drop)].index\n",
    "ls_type = [#'Sex',\n",
    "           'Organ',\n",
    "    #'Grade',\n",
    "    #'Stage',\n",
    "    #'Type',\n",
    "    #'Tissue ID.',\n",
    "    'Pathology',\n",
    "    'Age_Q']\n",
    "for s_type in ls_type:\n",
    "    print(s_type)\n",
    "    df_cluster[s_type] = df_cluster.scene.map(dict(zip(df_a.Scene,df_a.loc[:,s_type])))\n",
    "    confusion_matrix = pd.crosstab(df_cluster.loc[ls_index,'leiden'], df_cluster.loc[ls_index,s_type])\n",
    "    chi2, p, dof, exp = stats.chi2_contingency(confusion_matrix)\n",
    "    if p < 0.05:\n",
    "        print(f'{s_type} {p}')\n",
    "        #post hoc analysis\n",
    "        observed_vals = confusion_matrix\n",
    "        expected_vals = pd.DataFrame(exp,index=confusion_matrix.index,columns=confusion_matrix.columns)\n",
    "        col_sum = observed_vals.sum(axis=1)\n",
    "        row_sum = observed_vals.sum(axis=0)\n",
    "        leiden = list(confusion_matrix.index)\n",
    "        # Find all the unique attitudes towards flying for pairwise comparison\n",
    "        status = list(confusion_matrix.columns)\n",
    "        # Iterate through all combinations of parties and status\n",
    "        diff_vals = confusion_matrix - expected_vals\n",
    "        #annotated with pval\n",
    "        # cut out see cell below\n",
    "        cbar_kws = {'shrink':.3,'aspect':20,'anchor':(0.0, 1)}\n",
    "        fig, ax = plt.subplots(figsize = (5,5),dpi=300)\n",
    "        sns.heatmap(diff_vals.T,ax=ax,cmap='bwr',center=0,cbar_kws=cbar_kws)\n",
    "        ax.set_title(s_type)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'./Figures/confusion_matrix_difference_{s_sample}_{s_type}.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if s_type == 'Organ':\n",
    "    print('Organ')\n",
    "    matplotlib.rcParams.update({'font.size': 8})\n",
    "    #clustermap\n",
    "    fig = sns.clustermap(diff_vals.T,center=0,cmap='bwr',figsize = (5,6),  dendrogram_ratio=(.12, .06),\n",
    "               cbar_pos=None,#(0, .2, .03, .4)\n",
    "                        )\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/confusion_matrix_difference_{s_sample}_{s_type}_clustermap.png',dpi=300)\n",
    "if s_type == 'Pathology':\n",
    "    matplotlib.rcParams.update({'font.size': 8})\n",
    "    #clustermap\n",
    "    fig = sns.clustermap(diff_vals,center=0,cmap='bwr',figsize = (10,6),  dendrogram_ratio=(.08, .1),\n",
    "               cbar_pos=None,#(0, .2, .03, .4)\n",
    "                        )\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/confusion_matrix_difference_{s_sample}_{s_type}_clustermap.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        #annotated with pval\n",
    "'''\n",
    "        annot_vals = diff_vals.copy()\n",
    "        print(0.05/diff_vals.size)\n",
    "        for p in leiden:\n",
    "            for s in status:\n",
    "                observed = float(observed_vals.loc[p,s])\n",
    "                expected = float(expected_vals.loc[p,s])\n",
    "                col_total = float(col_sum[p])\n",
    "                row_total = float(row_sum[s])\n",
    "                expected_row_prop = expected/row_total\n",
    "                expected_col_prop = expected/col_total\n",
    "                std_resid = (observed - expected) / (math.sqrt(expected * (1-expected_row_prop) * (1-expected_col_prop)))\n",
    "                p_val = norm.sf(abs(std_resid))        \n",
    "                if p_val < 0.05/diff_vals.size:\n",
    "                    annot_vals.loc[p, s] =  \"*\"\n",
    "                else:\n",
    "                    annot_vals.loc[p, s] =  \"\"\n",
    "        #plot\n",
    "        annot_kws={'fontsize':7, \n",
    "           'fontstyle':'italic',  \n",
    "           'color':\"k\",\n",
    "           'alpha':0.6, \n",
    "        }\n",
    "        fig, ax = plt.subplots(figsize = (8,5),dpi=200)\n",
    "        sns.heatmap(diff_vals,ax=ax,annot=annot_vals,fmt='s',annot_kws=annot_kws,cmap='bwr',center=0)\n",
    "        '''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Percent H2O2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import intensity measurement and exposure time data\n",
    "ls_sample = ['4165NPanc','B1640','NPPan4165']\n",
    "df_data = pd.DataFrame()\n",
    "df_exp = pd.DataFrame()\n",
    "for s_sample in ls_sample:\n",
    "    print(s_sample)\n",
    "    df = pd.read_csv(f'{codedir}/Metadata/{s_sample}/{s_sample}_MeanIntensityMeasurement.csv',index_col=0)\n",
    "    df_data = df_data.append(df)\n",
    "    for s_tissue in df[df.tissue.str.contains(s_sample)].tissue:\n",
    "        df_exp_tissue = pd.read_csv(f'{codedir}/Metadata/{s_sample}/{s_tissue}_ExposureTimes.csv',index_col=0)\n",
    "        df_exp = df_exp.append(df_exp_tissue, sort='True')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#add experimental condition annotation\n",
    "d_condition = {'4165NPanc-73-Scene-001':'3%, 15min_stain',\n",
    " '4165NPanc-74-Scene-001':'3%, 15min_blank',\n",
    " '4165NPanc-75-Scene-001':'4.5%, 15min_stain',\n",
    " '4165NPanc-76-Scene-001':'4.5%, 15min_blank',\n",
    " '4165NPanc-77-Scene-001':'6%, 15min_stain',\n",
    " '4165NPanc-78-Scene-001':'6%, 15min_blank',\n",
    " 'B1640-T8-3-Scene-001':'3%, 30min_blank',\n",
    " 'B1640-T8-4-Scene-001':'3%, 60min_stain',\n",
    " 'B1640-T8-5-Scene-001':'3%, 30min_stain',\n",
    " 'B1640-T8-6-Scene-001':'3%, 30min+light_stain',\n",
    " 'NPPan4165-65-Scene-1': '3%, 30min+light_blank',\n",
    " 'NPPan4165-66-Scene-1': 'no-quench_blank',\n",
    "              }\n",
    "#add annotation\n",
    "for key, item in d_condition.items():\n",
    "    s_tissue = key.split('-Scene')[0]\n",
    "    s_type = item.split('_')[1]\n",
    "    s_cond = item.split('_')[0]\n",
    "    s_index = (df_data[df_data.scene==key]).index\n",
    "    df_data.loc[s_index,'condition'] = s_cond\n",
    "    df_data.loc[s_index,'slidetype'] = s_type\n",
    "    df_data.loc[s_index,'cond_type'] = item\n",
    "#normalize to 50 ms exposure time\n",
    "df_data['norm_result'] = (df_data.result/df_data.exposure)*50\n",
    "df_data = df_data.drop_duplicates(keep='first')\n",
    "d_color = {'c2':'AF488','c3':'AF555',\n",
    "           #'c4':'AF647','c5':'AF750'\n",
    "          }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for s_color, s_fluor in d_color.items():\n",
    "    fig,ax = plt.subplots(figsize=(5,3.5),dpi=200)\n",
    "    df_1 = df_data[(df_data.color==s_color) & (df_data.rounds=='R1') & (df_data.slidetype=='stain')].loc[:,['condition','norm_result']].set_index('condition')\n",
    "    df_2 = df_data[(df_data.color==s_color) & (df_data.rounds=='R2') & (df_data.slidetype=='stain')].loc[:,['condition','norm_result']].set_index('condition')\n",
    "    df_mean = df_data[(df_data.color==s_color) & (df_data.rounds!='R1') & (df_data.slidetype=='blank')].groupby('condition').mean().loc[:,'norm_result']\n",
    "    df_blank = (pd.DataFrame(df_mean)/df_1)\n",
    "    df_plot = (df_2/df_1).sort_values('norm_result',ascending=False)\n",
    "    #plot\n",
    "    df_plot.plot(legend=False,ax=ax)\n",
    "    if s_color == 'c2' or s_color == 'c3':\n",
    "        df_plot_b = df_blank.loc[df_blank.index.isin(df_plot.index)].fillna(df_blank.loc['3%, 30min','norm_result']).loc[df_plot.index]\n",
    "    else: \n",
    "        df_plot_b = df_blank.loc[df_blank.index.isin(df_plot.index)].loc[df_plot.index]\n",
    "    for tick in ax.get_xticklabels():\n",
    "        tick.set_rotation(45)\n",
    "    ax.set_ylabel('Signal Remaining')\n",
    "    ax.set_xlabel('H2O2 percent, time')\n",
    "    ax.set_title(f'{s_fluor} Signal Removal')\n",
    "    ax.set_ylim(0,1)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Quenching_Optimization_{s_color}.png')\n",
    "    df_plot_b.plot(label='blank',ax=ax,linestyle='--')\n",
    "    ax.legend(['stain','blank'])\n",
    "    for tick in ax.get_xticklabels():\n",
    "        tick.set_rotation(35)\n",
    "    fig.savefig(f'./Figures/Quenching_Optimization_{s_color}_plusAF.png')\n",
    " \n",
    "    #over af\n",
    "    fig,ax = plt.subplots(figsize=(5,3.5),dpi=200)\n",
    "    df_plot_both = (df_plot - df_plot_b).sort_values('norm_result',ascending=False)\n",
    "    df_plot_both.plot(ax=ax,legend=False)\n",
    "    for tick in ax.get_xticklabels():\n",
    "        tick.set_rotation(35)\n",
    "    ax.set_ylabel('Fraction Signal \\n Above Autoflourescence',fontsize=10)\n",
    "    ax.set_xlabel('H2O2 percent, time')\n",
    "    ax.set_title(f'{s_fluor} Signal Removal')\n",
    "    ax.set_ylim(-.1,.3)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./Figures/Quenching_Optimization_{s_color}_aboveAF.png')\n",
    "    #break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single Slide plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_tissue = {'4165NPanc':'Time'}\n",
    "d_color = {'c2':'AF488','c3':'AF555','c4':'AF647','c5':'AF750'}\n",
    "for s_tissue, s_test in d_tissue.items():\n",
    "    df_data_tissue = df_data[(df_data.tissue.str.find(s_tissue)==0) & (df_data.color!='c1')]\n",
    "    ls_color = sorted(set(df_data_tissue.color))\n",
    "    ls_cond = sorted(set(df_data_tissue.condition))\n",
    "    for s_cond in ls_cond:\n",
    "        fig,ax = plt.subplots(2,2, figsize=(7,5),squeeze=False,dpi=200)\n",
    "        ax=ax.ravel()\n",
    "        for idx,s_color in enumerate(ls_color):\n",
    "            ls_legend=[]\n",
    "            ls_tissue = sorted(set(df_data_tissue.tissue))\n",
    "            df_tissue_cond = df_data[(df_data.color==s_color)& (df_data.condition==s_cond)]\n",
    "            for s_type in ['stain','blank']:\n",
    "                df_tissue = df_tissue_cond[df_tissue_cond.slidetype==s_type]\n",
    "                df_tissue = df_tissue.sort_values('round_int')\n",
    "                s_label = sorted(set(df_tissue.loc[:,'cond_type']))[0]\n",
    "                df_tissue.norm_result.plot(ax=ax[idx])\n",
    "                ls_legend= ls_legend + [s_label]\n",
    "            ax[idx].set_title(d_color[s_color])\n",
    "            ax[idx].set_ylabel('Mean Intensity')\n",
    "            ax[idx].set_xlabel('Minutes of Quenching ')\n",
    "            ax[idx].set_xticks(list(df_tissue.round_int - 1))\n",
    "            ax[idx].set_xticklabels(list(((df_tissue.round_int - 1).astype(int)*15).astype(str)))\n",
    "            ax[idx].legend([item.replace('_',' ') for item in ls_legend],fontsize=10)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/vertical_{s_tissue}{s_cond}_quenching.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fraction quenched\n",
    "d_rename = {'CK7':'CK7-AF488','CoxIV':'CoxIV-AF555','LaminB1':'LamB1-AF647','Ecad':'Ecad-AF750'}\n",
    "for s_marker in ['CK7','CoxIV','LaminB1','Ecad']:\n",
    "    df_change = pd.DataFrame()\n",
    "    df_data[df_data.marker==s_marker]\n",
    "    for s_cond in ['3%','4.5%','6%']:\n",
    "        df = df_data[(df_data.marker==s_marker) & (df_data.condition==s_cond)]\n",
    "        df_change.loc[s_cond,'first'] = (df[df.rounds=='R1'].result[0] - df[df.rounds=='R2'].result[0])/(df[df.rounds=='R1'].result[0] - df[df.rounds=='R5'].result[0])\n",
    "        df_change.loc[s_cond,'second'] = (df[df.rounds=='R1'].result[0] - df[df.rounds=='R3'].result[0])/(df[df.rounds=='R1'].result[0] - df[df.rounds=='R5'].result[0])\n",
    "        df_change.loc[s_cond,'third'] = (df[df.rounds=='R1'].result[0] - df[df.rounds=='R4'].result[0])/(df[df.rounds=='R1'].result[0] - df[df.rounds=='R5'].result[0])\n",
    "    fig, ax = plt.subplots(figsize=(3.4,3),dpi=200)\n",
    "    df_change.plot(kind='bar',ax=ax,alpha=.8)\n",
    "    ax.legend(loc='lower right',title='quench no.',fontsize=10,title_fontsize=10)\n",
    "    ax.set_title(f'{d_rename[s_marker]}')\n",
    "    ax.set_ylabel('Fraction of Signal\\n Quenched')\n",
    "    ax.set_xlabel('H2O2 percentage')\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{codedir}/Figures/quenching_rate_by_H2O2_{s_marker}.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_quen = pd.DataFrame()\n",
    "#how much does AF reduce in 1st quench?\n",
    "for s_color in ['c2','c3','c4','c5']:\n",
    "    df_test = df_data[(df_data.slidetype=='blank') & (df_data.color==s_color) &(df_data.scene.str.contains('4165NPanc-')) & (df_data.rounds.isin(['R1','R2']))]\n",
    "    i_result = df_test.groupby('rounds').mean().result.loc['R2']/df_test.groupby('rounds').mean().result.loc['R1']\n",
    "    print(f'{s_color}: {1-i_result}')\n",
    "    ls_single = []\n",
    "    for s_tissue in ['4165NPanc-74', '4165NPanc-76', '4165NPanc-78']:\n",
    "        i_single = df_test[(df_test.tissue==s_tissue) & (df_test.rounds=='R2')].result[0]/(df_test[(df_test.tissue==s_tissue) & (df_test.rounds=='R1')].result[0])\n",
    "        ls_single.append(i_single)\n",
    "        df_quen.loc[f'{s_tissue}_{s_color}','result'] = 1 - i_single\n",
    "        df_quen.loc[f'{s_tissue}_{s_color}','channel'] = s_color\n",
    "    print(np.std(ls_single))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize= (3.2,3),dpi=200)\n",
    "sns.barplot(data=df_quen, x='channel', y='result', palette = 'muted',ci='sd',ax=ax)\n",
    "sns.stripplot(data=df_quen, x='channel', y='result', palette = 'dark',ax=ax)\n",
    "ax.set_ylabel('Percent')\n",
    "ax.set_title('Reduction in\\n Autofluorescence')\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/reduction_in_AF_quenching.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#how much does AF reduce in 1st quench?\n",
    "for s_color in ['c2','c3','c4','c5']:\n",
    "    df_test = df_data[(df_data.slidetype=='blank') & (df_data.color==s_color) &(df_data.scene.str.contains('4165NPanc-')) & (df_data.rounds.isin(['R1','R2']))]\n",
    "    i_result = df_test.groupby('rounds').std().result.loc['R2']/df_test.groupby('rounds').std().result.loc['R1']\n",
    "    print(f'{s_color}: {i_result}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just AF488\n",
    "d_tissue = {'4165NPanc':'Time'}\n",
    "d_color = {'c2':'AF488','c3':'AF555','c4':'AF647','c5':'AF750'}\n",
    "for s_tissue, s_test in d_tissue.items():\n",
    "    df_data_tissue = df_data[(df_data.tissue.str.find(s_tissue)==0) & (df_data.color!='c1')]\n",
    "    ls_color = sorted(set(df_data_tissue.color))\n",
    "    ls_cond = sorted(set(df_data_tissue.condition))\n",
    "    for s_cond in ls_cond:\n",
    "        fig,ax = plt.subplots(figsize=(4,3),squeeze=False,dpi=200)\n",
    "        ax=ax.ravel()\n",
    "        for idx,s_color in enumerate(ls_color):\n",
    "            ls_legend=[]\n",
    "            ls_tissue = sorted(set(df_data_tissue.tissue))\n",
    "            df_tissue_cond = df_data[(df_data.color==s_color)& (df_data.condition==s_cond)]\n",
    "            for s_type in ['stain','blank']:\n",
    "                df_tissue = df_tissue_cond[df_tissue_cond.slidetype==s_type]\n",
    "                df_tissue = df_tissue.sort_values('round_int')\n",
    "                s_label = sorted(set(df_tissue.loc[:,'cond_type']))[0]\n",
    "                df_tissue.norm_result.plot(ax=ax[idx])\n",
    "                ls_legend= ls_legend + [s_label]\n",
    "            ax[idx].set_title(d_color[s_color])\n",
    "            ax[idx].set_ylabel('Mean Intensity')\n",
    "            ax[idx].set_xlabel('Minutes of Quenching ')\n",
    "            ax[idx].set_xticks(list(df_tissue.round_int - 1))\n",
    "            ax[idx].set_xticklabels(list(((df_tissue.round_int - 1).astype(int)*15).astype(str)))\n",
    "            ax[idx].legend([item.replace('_',' ') for item in ls_legend],fontsize=12)\n",
    "            break\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/{s_tissue}/AF488_{s_tissue}{s_cond}_quenching.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Time and Addition of Incandescent Light"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_tissue = {'B1640':'Light'}\n",
    "d_color = {'c2':'AF488','c3':'AF555','c4':'AF647','c5':'AF750'}\n",
    "for s_tissue, s_test in d_tissue.items():\n",
    "    df_data_tissue = df_data[(df_data.tissue.str.contains(s_tissue)) & (df_data.color!='c1')]\n",
    "    ls_color = sorted(set(df_data_tissue.color))\n",
    "    ls_cond = ['30min', '60min', '30min+light']\n",
    "    for s_color in ls_color:\n",
    "        fig,ax = plt.subplots(2,2, figsize=(7,5),squeeze=False,dpi=200)\n",
    "        ax=ax.ravel()\n",
    "        for idx,s_cond in enumerate(ls_cond):\n",
    "            ls_legend=[]\n",
    "            ls_tissue = sorted(set(df_data_tissue.tissue))\n",
    "            df_tissue_cond = df_data[(df_data.color==s_color)& (df_data.condition==s_cond)]\n",
    "            for s_type in ['stain','blank']:\n",
    "                if s_type == 'blank':\n",
    "                    df_tissue = df_data[(df_data.cond_type=='30min_blank') & (df_data.color==s_color)]\n",
    "                else:\n",
    "                    df_tissue = df_tissue_cond[df_tissue_cond.slidetype==s_type]\n",
    "                df_tissue = df_tissue.sort_values('round_int')\n",
    "                s_label = sorted(set(df_tissue.loc[:,'cond_type']))[0]\n",
    "                df_tissue.norm_result.plot(ax=ax[idx])\n",
    "                ls_legend= ls_legend + [s_label]\n",
    "            ax[idx].set_title(d_color[s_color])\n",
    "            ax[idx].set_ylabel('Mean Intensity')\n",
    "            ax[idx].set_xlabel('Minutes of Quenching')\n",
    "            ax[idx].set_xticks(list(df_tissue.round_int - 1))\n",
    "            fontdict = {'fontsize': 10}\n",
    "            ax[idx].set_xticklabels(list(((df_tissue.round_int - 1).astype(int)*30).astype(str)),fontdict=fontdict)\n",
    "            if idx==1:\n",
    "                ax[idx].set_xticklabels(list(((df_tissue.round_int - 1).astype(int)*60).astype(str)),fontdict=fontdict)\n",
    "            ax[idx].legend([item.replace('_',' ') for item in ls_legend],fontsize=10)\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/{s_tissue}/vertical_{s_tissue}{s_color}_quenching.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#just AF488\n",
    "d_tissue = {'B1640':'Light'}\n",
    "d_color = {'c2':'AF488','c3':'AF555','c4':'AF647','c5':'AF750'}\n",
    "for s_tissue, s_test in d_tissue.items():\n",
    "    df_data_tissue = df_data[(df_data.tissue.str.contains(s_tissue)) & (df_data.color!='c1')]\n",
    "    ls_color = sorted(set(df_data_tissue.color))\n",
    "    ls_cond = ['30min', '60min', '30min+light']\n",
    "    for s_cond in ls_cond:\n",
    "        fig,ax = plt.subplots(figsize=(4,3),squeeze=False,dpi=200)\n",
    "        ax=ax.ravel()\n",
    "        for idx,s_color in enumerate(ls_color):\n",
    "            ls_legend=[]\n",
    "            ls_tissue = sorted(set(df_data_tissue.tissue))\n",
    "            df_tissue_cond = df_data[(df_data.color==s_color)& (df_data.condition==s_cond)]\n",
    "            for s_type in ['stain','blank']:\n",
    "                if s_type == 'blank':\n",
    "                    df_tissue = df_data[(df_data.cond_type=='30min_blank') & (df_data.color==s_color)]\n",
    "                else:\n",
    "                    df_tissue = df_tissue_cond[df_tissue_cond.slidetype==s_type]\n",
    "                df_tissue = df_tissue.sort_values('round_int')\n",
    "                s_label = sorted(set(df_tissue.loc[:,'cond_type']))[0]\n",
    "                df_tissue.norm_result.plot(ax=ax[idx])\n",
    "                ls_legend= ls_legend + [s_label]\n",
    "            ax[idx].set_title(d_color[s_color])\n",
    "            ax[idx].set_ylabel('Mean Intensity')\n",
    "            ax[idx].set_xlabel('Minutes of Quenching')\n",
    "            ax[idx].set_xticks(list(df_tissue.round_int - 1))\n",
    "            fontdict = {'fontsize': 10}\n",
    "            ax[idx].set_xticklabels(list(((df_tissue.round_int - 1).astype(int)*30).astype(str)),fontdict=fontdict)\n",
    "            if idx==1:\n",
    "                ax[idx].set_xticklabels(list(((df_tissue.round_int - 1).astype(int)*60).astype(str)),fontdict=fontdict)\n",
    "            ax[idx].legend([item.replace('_',' ') for item in ls_legend],fontsize=12)\n",
    "            break\n",
    "        plt.tight_layout()\n",
    "        fig.savefig(f'{codedir}/Figures/{s_tissue}/AF488_{s_tissue}{s_color}{s_cond}_quenching.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Blank Slide Quenched and Imaged vs. Imaged Only\n",
    "\n",
    "Does the reducion in autofluorescnce come from the quenching or imaging process (i.e. is it chemical or photo bleaching of autofluorescnce) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_tissue = {'NPPan4165':'No Quench Ctrl'}\n",
    "d_color = {'c2':'AF488','c3':'AF555','c4':'AF647','c5':'AF750'}\n",
    "d_label = {'no-quench_blank':'imaging only','30min+light_blank':'quench + image'}\n",
    "for s_tissue, s_test in d_tissue.items():\n",
    "    df_data_tissue = df_data[(df_data.tissue.str.contains(s_tissue)) & (df_data.color!='c1')]\n",
    "    ls_color = sorted(set(df_data_tissue.color))\n",
    "    fig,ax = plt.subplots(2,2, figsize=(7,5.5),squeeze=False,dpi=200)\n",
    "    ax=ax.ravel()\n",
    "    for idx,s_color in enumerate(ls_color):\n",
    "        ls_legend=[]\n",
    "        ls_tissue = sorted(set(df_data_tissue.tissue))\n",
    "        for s_tissue in ls_tissue:\n",
    "            df_tissue = df_data[(df_data.tissue==s_tissue) & (df_data.color==s_color)]\n",
    "            df_tissue = df_tissue.sort_values('round_int')\n",
    "            s_label = d_label[sorted(set(df_tissue.loc[:,'cond_type']))[0]]\n",
    "            df_tissue.norm_result.plot(ax=ax[idx])\n",
    "            ls_legend= ls_legend + [s_label]\n",
    "        ax[idx].set_title(d_color[s_color])\n",
    "        ax[idx].set_ylabel('Mean Intensity')\n",
    "        ax[idx].set_xlabel('Minutes of Quenching')\n",
    "        ax[idx].set_xticks(list(df_tissue.round_int - 1))\n",
    "        ax[idx].set_xticklabels(list(((df_tissue.round_int - 1)*30).astype(int).astype(str)))\n",
    "        ax[idx].legend([item.replace('_',' ') for item in ls_legend],fontsize=8)\n",
    "        ax[idx].set_xlim(left=0)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{codedir}/Figures/NPPan4165/vertical_{s_tissue}_imagingonly.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#just AF488\n",
    "d_tissue = {'NPPan4165':'No Quench Ctrl'}\n",
    "d_color = {'c2':'AF488','c3':'AF555','c4':'AF647','c5':'AF750'}\n",
    "d_label = {'no-quench_blank':'imaging only','30min+light_blank':'quench + image'}\n",
    "for s_tissue, s_test in d_tissue.items():\n",
    "    df_data_tissue = df_data[(df_data.tissue.str.contains(s_tissue)) & (df_data.color!='c1')]\n",
    "    ls_color = sorted(set(df_data_tissue.color))\n",
    "    fig,ax = plt.subplots(figsize=(4,3),squeeze=False,dpi=200)\n",
    "    ax=ax.ravel()\n",
    "    for idx,s_color in enumerate(ls_color):\n",
    "        ls_legend=[]\n",
    "        ls_tissue = sorted(set(df_data_tissue.tissue))\n",
    "        for s_tissue in ls_tissue:\n",
    "            df_tissue = df_data[(df_data.tissue==s_tissue) & (df_data.color==s_color)]\n",
    "            df_tissue = df_tissue.sort_values('round_int')\n",
    "            s_label = d_label[sorted(set(df_tissue.loc[:,'cond_type']))[0]]\n",
    "            df_tissue.norm_result.plot(ax=ax[idx])\n",
    "            ls_legend= ls_legend + [s_label]\n",
    "        ax[idx].set_title(d_color[s_color])\n",
    "        ax[idx].set_ylabel('Mean Intensity')\n",
    "        ax[idx].set_xlabel('Minutes of Quenching')\n",
    "        ax[idx].set_xticks(list(df_tissue.round_int - 1))\n",
    "        ax[idx].set_xticklabels(list(((df_tissue.round_int - 1)*30).astype(int).astype(str)))\n",
    "        ax[idx].legend([item.replace('_',' ') for item in ls_legend],fontsize=12)\n",
    "        ax[idx].set_xlim(left=0)\n",
    "        break\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{codedir}/Figures/NPPan4165/AF488_{s_tissue}_imagingonly.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Temperature of Buffer under Various Conditions\n",
    "\n",
    "Our intepretation is that the heat comes from the light, but also from the exotermic reaction of H2O2 oxidation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_temp = pd.read_csv(f'{codedir}/Data/TemperatureExperiments.csv',index_col = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 12})\n",
    "fig, ax = plt.subplots(dpi=200)\n",
    "df_temp.iloc[:7,[0,2,4]].plot(kind='line', ax=ax)\n",
    "ax.set_ylabel('Temp (C)')\n",
    "ax.set_xlabel('Minutes')\n",
    "ax.legend(labels = [item.split('_')[0] for item in ax.get_legend_handles_labels()[1]],fontsize=10)\n",
    "ax.set_title('Effect of Light on Solution Temperature')\n",
    "fig.savefig(f'{codedir}/Figures/EffectofLightonQuenchTemperature.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_temp.iloc[:7,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 12})\n",
    "fig, ax = plt.subplots(dpi=200)\n",
    "df_temp.iloc[:7,[0,1,2,3]].plot(kind='line', ax=ax)\n",
    "ax.set_ylabel('Temp (C)')\n",
    "ax.set_xlabel('Minutes')\n",
    "ax.legend(labels = [item.split('_')[0] + item.split('_')[2] for item in ax.get_legend_handles_labels()[1]],fontsize=10)\n",
    "ax.set_title('Effect of H2O2 on Solution Temperature')\n",
    "fig.savefig(f'{codedir}/Figures/EffectofH2O2onQuenchTemperature.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Temperature result interpretation\n",
    "\n",
    "The incandescent light physically heats the solution, but also increases the rate of H2O2 oxidation and exothermic heat production. This effect is not seen with LED light.\n",
    "\n",
    "We also observed the most bubbles in the Incandescent/H2O2 solution, suggesting the highest level of oxidation."
   ]
  },
  {
   "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
}