Normalization_testing_tissue.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Batch Normalization on tissues\n",
    "\n",
    "\n",
    "**Question:** How does combat batch normalization algorithms perform on adjacent section tissues stained with the same panel, when training with unmatched cores?\n",
    "\n",
    "**Samples:** \n",
    "- TMA: Purchased from biomax: https://www.biomax.us/tissue-arrays/Breast/BR1506\n",
    "- Adjacent Sections: BM-Her2N75-15, BM-Her2N75-17, BM-Her2N75-18 (section 16 skipped.)\n",
    "- Scenes: (i.e. TMA cores) \n",
    "  - 17: ER+/HER2+, immune rich. \n",
    "  - 49: ER+/HER2+\n",
    "  - 59: HER2+ immune rich\n",
    "\n",
    "**Method**: We performed combat normalization using unlike tissues as training set (different cores in training set and all cores in testing set). We visualized the resulting histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load libraries\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import copy\n",
    "import seaborn as sns\n",
    "import scipy\n",
    "import scanpy as sc\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.preprocessing import scale, minmax_scale, StandardScaler\n",
    "from sklearn.metrics import silhouette_score\n",
    "import matplotlib as mpl\n",
    "import util\n",
    "mpl.rc('figure', max_open_warning = 0)\n",
    "#os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cmIF')\n",
    "#from mplex_image import visualize as viz, process, preprocess, normalize\n",
    "np.random.seed(1202)\n",
    "np.random.seed(1211)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#change to correct directory\n",
    "os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cycIF_ValidationStudies/cycIF_Validation')\n",
    "codedir=os.getcwd()\n",
    "rootdir = f'{codedir}/Data/'\n",
    "datadir = f'{codedir}/Data/filtered_data'\n",
    "os.chdir(datadir)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "df=pd.read_csv(f'20201229_BM-Her2N75-15-17-18_MeanIntensity.csv',index_col=0)\n",
    "df.rename({'slide':'batch'},axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.groupby('slide_scene').mean().loc[:,['CK7_Ring','CK14_Ring','CK5_Ring','CK19_Ring','Ecad_Ring','HER2_Ring']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#all cores\n",
    "data = df.loc[:,df.dtypes=='float64'].T\n",
    "batch = df.batch\n",
    "gamma_star, delta_star, stand_mean, var_pooled = util.combat_fit(data, batch)\n",
    "#transform\n",
    "bayesdata = util.combat_transform(data,batch,gamma_star, delta_star, stand_mean, var_pooled)\n",
    "df_norm=bayesdata.T\n",
    "s_train='all'\n",
    "s_tissue = 'all'\n",
    "df_norm['batch'] = [item.split('_')[0] for item in df_norm.index]\n",
    "d_fig = util.plot_histograms(df_norm,df,s_train,s_tissue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_marker, fig in d_fig.items():\n",
    "    fig.savefig(f'{rootdir}/20201228/Different_Scaling_combat_training_{s_train}_{s_marker}_{s_tissue}.png')\n",
    "    fig\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training combat with controls\n",
    "\n",
    "How sensitive is the combat algorithm to the inputs that determine the parameters?\n",
    "\n",
    "- Can different tissues in each batch be used to fit combat?\n",
    "- Can a limited set of tissues be used to fit combat?\n",
    "- Or, should the set of tissues used to fit combat be very similar to those it is applied to?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lls_batch = [['BM-Her2N75-15', 'BM-Her2N75-17', 'BM-Her2N75-18'],\n",
    "            [ 'BM-Her2N75-17', 'BM-Her2N75-18','BM-Her2N75-15'],\n",
    "            [ 'BM-Her2N75-18','BM-Her2N75-15', 'BM-Her2N75-17']]\n",
    "ls_scene = sorted(set(df.scene))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#normalize with different training sets\n",
    "for idxx, ls_batch in enumerate(lls_batch):\n",
    "    data = pd.DataFrame()\n",
    "    for idx, s_batch in enumerate(ls_batch):\n",
    "        s_scene = ['scene017', 'scene049', 'scene059'][idx]\n",
    "        data = data.append(df.loc[((df.scene==s_scene)&(df.batch==s_batch)),:])   \n",
    "    #fit training set\n",
    "    gamma_star, delta_star, stand_mean, var_pooled = util.combat_fit(data.loc[:,data.dtypes=='float64'].T, data.batch)\n",
    "    #transform full data set\n",
    "    bayesdata = util.combat_transform(df.loc[:,df.dtypes=='float64'].T,df.batch,gamma_star, delta_star,stand_mean, var_pooled)\n",
    "    df_norm=bayesdata.T\n",
    "    s_train = \"_\".join([(item + ls_scene[idx]).split('-')[-1].replace('scene','s') for idx, item in enumerate(ls_batch)])\n",
    "    s_tissue = 'diff'\n",
    "    d_fig = util.plot_histograms(df_norm,data,s_train,s_tissue)\n",
    "    #break\n",
    "    df_norm.sample(5400,random_state=3).to_csv(f'20210301-{idxx}_BM-Her2N75_SampledMeanIntensity_diff_train.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for s_marker, fig in d_fig.items():\n",
    "    fig.savefig(f'{rootdir}/20201228/Different_Scaling_combat_training_{s_train}_{s_marker}_{s_tissue}.png')\n",
    "    fig\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#normalize with same training sets\n",
    "for idx, ls_batch in enumerate(lls_batch):\n",
    "    s_scene = ['scene017', 'scene049', 'scene059'][idx]\n",
    "    data = pd.DataFrame()\n",
    "    for s_batch in ls_batch:   \n",
    "        data = data.append(df.loc[((df.scene==s_scene)&(df.batch==s_batch)),:])   \n",
    "    #fit training set\n",
    "    gamma_star, delta_star, stand_mean, var_pooled = util.combat_fit(data.loc[:,data.dtypes=='float64'].T, data.batch)\n",
    "    #transform full data set\n",
    "    bayesdata = util.combat_transform(df.loc[:,df.dtypes=='float64'].T,df.batch,gamma_star, delta_star, stand_mean, var_pooled)\n",
    "    df_norm=bayesdata.T\n",
    "    s_train = \"_\".join([(item + s_scene).split('-')[-1].replace('scene','s') for idx, item in enumerate(ls_batch)])\n",
    "    s_tissue = 'same'\n",
    "    d_fig = util.plot_histograms(df_norm,data,s_train,s_tissue)\n",
    "    break\n",
    "    df_norm.sample(5400,random_state=3).to_csv(f'20210301-{idx}_BM-Her2N75_SampledMeanIntensity_same_train.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_marker, fig in d_fig.items():\n",
    "    fig.savefig(f'{rootdir}/20201228/Different_Scaling_combat_training_{s_train}_{s_marker}_{s_tissue}.png')\n",
    "    fig\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note: \n",
    "20210301-0_BM-Her2N75-15-17-18_kbet_same_train.csv gave the worst kbet results\n",
    "this is the ER+/HER2+, immune rich core\n",
    "markers failed on were CK14,  CK5, CK7, and Ecad; HER2 and CK19 to a lesser degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#normalize with sampled training sets\n",
    "ls_date = ['20201207','20201208','20201209']\n",
    "for idx, s_date in enumerate(ls_date):\n",
    "    data = pd.read_csv(f'{s_date}_BM-Her2N75-15-17-18_SampledMeanIntensity_raw.csv',index_col=0)\n",
    "    data['batch'] = [item.split('_')[0] for item in data.index]\n",
    "    #fit training set\n",
    "    gamma_star, delta_star, stand_mean, var_pooled = util.combat_fit(data.loc[:,data.dtypes=='float64'].T, data.batch)\n",
    "    #transform full data set\n",
    "    bayesdata = util.combat_transform(df.loc[:,df.dtypes=='float64'].T,df.batch,gamma_star, delta_star, stand_mean, var_pooled)\n",
    "    df_norm=bayesdata.T\n",
    "    s_train = s_date\n",
    "    s_tissue = 'sampled'\n",
    "    d_fig =util.plot_histograms(df_norm,data,s_train,s_tissue)\n",
    "    break\n",
    "    df_norm.sample(5400,random_state=3).to_csv(f'20210301-{idx}_BM-Her2N75_SampledMeanIntensity_sampled_train.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for s_marker, fig in d_fig.items():\n",
    "    fig.savefig(f'{rootdir}/20201228/Different_Scaling_combat_training_{s_train}_{s_marker}_{s_tissue}.png')\n",
    "    fig\n",
    "    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RESTORE methods: Results\n",
    "\n",
    "What is the best practice for restore normalization?\n",
    "\n",
    "- Global or local thresholds? - **Local slightly better but neither very good**\n",
    "- Raw data or arcsinh transformation? - **No difference**\n",
    "- Divide by threshold or scale values above threshold? - **Scale gives good normaliztion, but not as good as combat**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_file = pd.DataFrame(index=os.listdir())\n",
    "df_file = df_file[df_file.index.str.contains('BM-Her2N75-15-17-18_kbet')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir(filterdir)\n",
    "df_file = pd.DataFrame(index=os.listdir())\n",
    "df_file = df_file[df_file.index.str.contains('BM-Her2N75-15-17-18_kbet_')]\n",
    "#ls_train = ['diff_train','raw', 'same_train', 'sampled_train','raw_combat']\n",
    "ls_method = ['raw', 'raw_combat', 'raw_regress_out', 'restore_local', 'restore_div','restore_div_arcsinh', 'restore_scale','restore_scale_combat']#,'raw_restore_combat' 'raw_restore','raw_restore_regress_out',\n",
    "\n",
    "#add mean kbet\n",
    "for s_file in df_file.index:\n",
    "    df = pd.read_csv(s_file,index_col=0)\n",
    "    df_file.loc[s_file,'mean_kbet'] = df.loc['mean','kBET.observed']\n",
    "df_file['norm'] = [item.split('kbet_')[1].split('.csv')[0] for item in df_file.index]\n",
    "df_file = df_file[df_file.norm.isin(ls_method)]\n",
    "ls_index= df_file.groupby('norm').mean_kbet.mean().sort_values().index\n",
    "df_file.groupby('norm').mean_kbet.mean().sort_values()\n",
    "df_file['Norm'] = pd.Categorical(\n",
    "    df_file['norm'], \n",
    "    categories=ls_index.tolist(), \n",
    "    ordered=True\n",
    ")\n",
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(5,3),dpi=200)\n",
    "df_plot = df_file.sort_values('Norm')\n",
    "sns.lineplot(data=df_plot,x='norm',y='mean_kbet',ax=ax,err_style='bars')\n",
    "labels = [item.replace('raw_','').replace('_','\\n').replace('div','global') for item in ls_index.tolist()]\n",
    "ax.set_xticks(range(len(df_file.groupby('norm').mean_kbet)))\n",
    "ax.set_xticklabels(labels,rotation=0)\n",
    "ax.set_ylabel('Rejection Rate')\n",
    "ax.set_xlabel('Normalization')\n",
    "ax.set_title('kBET Evaluation of Batch Correction II')\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "fig.savefig(f'{rootdir}/20201228/BatchEffectII.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot.groupby('Norm').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot.groupby('Norm').std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_select = df_plot[df_plot.Norm.isin(['raw_combat','restore_scale','raw_regress_out','raw'])]\n",
    "fig, ax = plt.subplots(figsize=(4.1,3),dpi=300)\n",
    "df_plot = df_select.sort_values('Norm')\n",
    "sns.lineplot(data=df_plot,x='norm',y='mean_kbet',ax=ax,err_style='bars')\n",
    "labels = ['combat','RESTORE','regress\\n out','raw']\n",
    "ax.set_xticks(range(len(df_select.groupby('norm').mean_kbet)))\n",
    "ax.set_xticklabels(labels,rotation=0,fontsize=13)\n",
    "ax.set_ylabel('Rejection Rate',fontsize=14)\n",
    "ax.set_xlabel('Normalization',fontsize=14)\n",
    "ax.set_ylim(.4,1)\n",
    "ax.set_title('kBET Evaluation of\\n Batch Correction', fontsize=18)\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "fig.savefig(f'{codedir}/Figures/BatchEffectII_select.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"whitegrid\")\n",
    "df_pearson = pd.read_csv(f'Pearson_correlation_tissue_0.6.csv',index_col=0)\n",
    "ls_order = df_pearson.mean(axis=1).sort_values(ascending=False).index.tolist()\n",
    "#ls_order = ['raw_combat','log2_combat','restore_scale','raw_regress_out','raw','raw_robust','raw_standard']\n",
    "fig,ax = plt.subplots(figsize=(3.5,2.6),dpi=300)\n",
    "sns.boxplot(data=df_pearson.loc[ls_order].T,ax=ax,orient='h',showfliers=False,palette='muted')\n",
    "sns.stripplot(data=df_pearson.loc[ls_order].T,ax=ax,orient='h',palette='dark')\n",
    "#ax.set_title(f'Cluster Correlation \\n (resolution {resolution})',fontsize=16)\n",
    "ax.set_title(f'Cluster Correlation',fontsize=16)\n",
    "ax.set_xlabel('Pearson Correlation',fontsize=14)\n",
    "ax.yaxis.set_label_position(\"right\")\n",
    "ax.set_yticklabels(['combat','RESTORE','regress out','raw'],fontsize=14)\n",
    "ax.yaxis.tick_right()\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{rootdir}filtered_data/figures/PearsonCorrelation_tissue_0.6.png',dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training combat with controls: Results\n",
    "\n",
    "How sensitive is the combat algorithm to the inputs that determine the parameters?\n",
    "\n",
    "- Can different tissues in each batch be used to fit combat? - **No, kbet rejection rate is near 1**\n",
    "- Can a limited set of tissues be used to fit combat? - **Varies, kbet rejection from 0.85 to 1, depending on training tissue**\n",
    "- Or, should the set of tissues used to fit combat be very similar to those it is applied to? - **Yes, lowest kbet rejection rate**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#add mean kbet\n",
    "df_select = df_file[df_file.norm.isin(['raw_combat','raw_standard','restore_scale','raw','raw_regress_out','raw_robust','log2_combat'])] #'raw_restore_combat'\n",
    "sns.set_style(\"white\")\n",
    "df_file = pd.DataFrame(index=os.listdir())\n",
    "df_file = df_file[df_file.index.str.contains('BM-Her2N75-15-17-18_kbet_')]\n",
    "ls_train = ['diff_train','raw', 'same_train', 'sampled_train','raw_combat']\n",
    "for s_file in df_file.index:\n",
    "    df = pd.read_csv(s_file,index_col=0)\n",
    "    df_file.loc[s_file,'mean_kbet'] = df.loc['mean','kBET.observed']\n",
    "df_file['norm'] = [item.split('kbet_')[1].split('.csv')[0] for item in df_file.index]\n",
    "df_file = df_file[df_file.norm.isin(ls_train)]\n",
    "ls_index= df_file.groupby('norm').mean_kbet.mean().sort_values().index\n",
    "df_file.groupby('norm').mean_kbet.mean().sort_values()\n",
    "df_file['Norm'] = pd.Categorical(\n",
    "    df_file['norm'], \n",
    "    categories=ls_index.tolist(), \n",
    "    ordered=True\n",
    ")\n",
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(3.8,3),dpi=300)\n",
    "df_plot = df_file.sort_values('Norm')\n",
    "sns.lineplot(data=df_plot,x='norm',y='mean_kbet',ax=ax,err_style='bars')\n",
    "labels = [item.replace('raw_','baseline_').replace('_','\\n').replace('train','training').replace('diff','different') for item in ls_index.tolist()]\n",
    "ax.set_xticks(range(len(df_file.groupby('norm').mean_kbet)))\n",
    "ax.set_xticklabels(labels,rotation=90,fontsize=13)\n",
    "ax.set_ylabel('Rejection Rate',fontsize=14)\n",
    "ax.set_xlabel('Normalization',fontsize=14)\n",
    "ax.set_ylim(.4,1)\n",
    "ax.set_title('kBET Evaluation of\\n Combat Parameterization',fontsize=18,x=.43)\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "fig.savefig(f'{codedir}/Figures/BatchEffect_CombatParam.png',dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot.groupby('Norm').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot.groupby('Norm').std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_file.sort_values('Norm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_index = [ '20201208_BM-Her2N75-15-17-18_SampledMeanIntensity_raw.csv',\n",
    " '20210301-2_BM-Her2N75_SampledMeanIntensity_same_train.csv',\n",
    " '20210301-2_BM-Her2N75_SampledMeanIntensity_diff_train.csv',\n",
    " '20210301-0_BM-Her2N75_SampledMeanIntensity_sampled_train.csv',\n",
    " '20201209_BM-Her2N75-15-17-18_SampledMeanIntensity_raw_combat.csv'\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "marker_genes = ['CD44', 'CD45', 'CD4', 'CD68', 'CK14', 'CK19', 'CK5', 'CK7', 'ER',\n",
    "       'Ecad', 'Ki67', 'PD1', 'Vim', 'aSMA', 'pHH3']\n",
    "d_je_tma = {'scene01':'tonsil1','scene02':'HCC1143', 'scene03':'HCC3153', 'scene04':'NBreast',\n",
    "    'scene05':'T47D','scene06':'T47D','scene07':'tonsil2','scene08':'BT474','scene09':'BT474','scene10':'AU565',\n",
    "    'scene11':'AU565','scene12':'MDAMB-436','scene13':'MDAMB-436' }\n",
    "d_je_tma = {'scene017': 'ER-HER2-imm',\n",
    "  'scene049': 'ER-HER2',\n",
    "  'scene059': 'HER2-imm'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalize \n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.metrics import adjusted_rand_score\n",
    "%matplotlib inline\n",
    "#sc.set_figure_params(scanpy=True, fontsize=14)\n",
    "df_pearson = pd.DataFrame()\n",
    "for s_index in ls_index:\n",
    "    print(s_index)\n",
    "    df = pd.read_csv(s_index,index_col=0)\n",
    "    df.columns = [item.split('_')[0] for item in df.columns]\n",
    "    marker_genes = df.columns[df.dtypes=='float64'].tolist()\n",
    "    adata = sc.AnnData(df.loc[:,df.dtypes=='float64'])\n",
    "    adata.obs['batch'] = [item.split('_scene')[0].split('BM-Her2N75-')[1] for item in adata.obs.index]\n",
    "    adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]\n",
    "    s_norm = ''#s_index.split('SampledMeanIntensity_')[1].split('.')[0]\n",
    "    adata.raw = adata\n",
    "    #reduce dimensionality\n",
    "    sc.tl.pca(adata, svd_solver='auto')\n",
    "    #save normalized data\n",
    "    df = pd.DataFrame(data=adata.X,index=adata.obs.index,columns=adata.var.index)\n",
    "    df['batch'] = [item.split('_scene')[0] for item in df.index]\n",
    "    #df.to_csv(f'{s_index.replace(\".csv\",f\"_{s_norm}.csv\")}')\n",
    "    # calculate neighbors     \n",
    "    sc.pp.neighbors(adata, n_neighbors=10, n_pcs=15)\n",
    "    sc.tl.umap(adata)\n",
    "    #umap plot\n",
    "    s_type = s_index.split('SampledMeanIntensity_')[1].split('.')[0]\n",
    "    figname = f\"UmapBatch_param_{s_type}_{s_norm}.png\"\n",
    "    fig,ax = plt.subplots(figsize=(5,5), dpi=200)\n",
    "    sc.pl.umap(adata, color='batch',title=f\"{s_type.replace('_',' ')} {s_norm}\",wspace=.25,ax=ax,save=figname)\n",
    "    fig,ax = plt.subplots(figsize=(5,5), dpi=200)\n",
    "    figname = f'UmapScene_param_{s_type}_{s_norm}.png'\n",
    "    fig = sc.pl.umap(adata, color='scene',save=figname,ax=ax,title=f\"{s_type.replace('_',' ')} {s_norm}\")\n",
    "    X_pca = adata.obsm['X_pca'] \n",
    "    #leiden\n",
    "    #resolution=0.25\n",
    "    resolution=0.4\n",
    "    #resolution=0.6\n",
    "    sc.tl.leiden(adata,resolution=resolution)\n",
    "    fig,ax = plt.subplots(figsize=(4.5,5),dpi=200)\n",
    "    figname=f'leiden_param_{s_type}_{s_norm}{resolution}.png'\n",
    "    sc.pl.umap(adata, color='leiden',ax=ax,save=figname)\n",
    "    fig,ax = plt.subplots(figsize=(6,3.7), dpi=200)\n",
    "    figname=f'Matrixplot_param_leiden_{s_type}_{s_norm}{resolution}.png'\n",
    "    sc.pl.matrixplot(adata, var_names=marker_genes, groupby=f'leiden',title=s_type.replace('_',' '),\n",
    "                     dendrogram=True,ax=ax,save=figname,standard_scale='var',colorbar_title='Relative\\nintensity')\n",
    "    df = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "    df['leiden'] = adata.obs['leiden']\n",
    "    #stacked bar\n",
    "    s_trans = s_index.split('Intensity_')[1].split('.')[0]\n",
    "    df['slide'] = [item.split('_')[0] for item in df.index]\n",
    "    df['scene'] = [d_je_tma[item.split('_')[1]] for item in df.index]\n",
    "    df['slide_scene'] = df.slide + '_' + df.scene\n",
    "    df_prop = (df.groupby([f'leiden','slide_scene']).CD4.count())/(df.groupby(['slide_scene']).CD4.count())\n",
    "    df_prop = df_prop.unstack().fillna(value=0).T\n",
    "    #barplot\n",
    "    fig,ax=plt.subplots(figsize=(5,3.7), dpi=200)\n",
    "    df_prop.columns = df_prop.columns.add_categories(['slide','scene'])\n",
    "    df_prop.index = [item.replace('BM-Her2N75-','') for item in df_prop.index]\n",
    "    df_prop['slide'] =[item.split('_')[0] for item in df_prop.index]\n",
    "    df_prop['scene'] =[item.split('_')[1] for item in df_prop.index]\n",
    "    df_prop.sort_values(['scene','slide']).plot(kind='bar',stacked=True,ax=ax,legend=True,cmap='tab20',width=.8)\n",
    "    ax.legend(bbox_to_anchor=(1.02, 1.2), ncol=1,fontsize=12)\n",
    "    ax.set_ylabel('Fraction Positive')\n",
    "    labels = ax.get_xticklabels()\n",
    "    ax.set_xticklabels(labels, fontsize=14)\n",
    "    ax.set_title(f\"{s_trans.replace('_',' ')} {s_norm}\")\n",
    "    ax.grid(False)\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'./figures/StackedBar_param_{s_trans}_{s_norm}{resolution}_Leiden.png')\n",
    "    #pearson\n",
    "    ls_core = sorted(set(df_prop.scene))\n",
    "    se_all = pd.Series(dtype='float64',name=s_trans)\n",
    "    for s_core in ls_core:\n",
    "        test = df_prop.loc[df_prop.scene==s_core,df_prop.dtypes=='float64']\n",
    "        se_test = pd.Series([scipy.stats.pearsonr(test.iloc[0],test.iloc[1])[0],scipy.stats.pearsonr(test.iloc[1],test.iloc[2])[0],scipy.stats.pearsonr(test.iloc[2],test.iloc[0])[0]])\n",
    "        se_all = se_all.append(se_test)\n",
    "    print(se_all.mean())\n",
    "    se_all.name = s_trans\n",
    "    se_all.index = (range(len(se_all)))\n",
    "    df_pearson = df_pearson.append(se_all)\n",
    "    #break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_pearson.to_csv(f'20210301_Pearson_correlation_param{resolution}.csv')\n",
    "rootdir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 0.6\n",
    "#df_pearson.to_csv(f'Pearson_correlation_param{resolution}.csv')\n",
    "df_pearson = pd.read_csv(f'20210301_Pearson_correlation_param{resolution}.csv',index_col=0)\n",
    "ls_order = df_pearson.mean(axis=1).sort_values(ascending=False).index.tolist()\n",
    "%matplotlib inline\n",
    "s_date = '20201228'\n",
    "fig,ax = plt.subplots(figsize=(4,4),dpi=200)\n",
    "sns.boxplot(data=df_pearson.loc[ls_order].T,ax=ax,orient='h',showfliers=False,palette='muted')\n",
    "sns.stripplot(data=df_pearson.loc[ls_order].T,ax=ax,orient='h',palette='dark')\n",
    "ax.set_title(f'Cluster Correlation \\n (resolution {resolution})',fontsize=16)\n",
    "#ax.set_title(f'Cluster Correlation',fontsize=16)\n",
    "ax.set_xlabel('Pearson Correlation',fontsize=14)\n",
    "if resolution == 0.6:\n",
    "    #pass\n",
    "    ax.set_yticklabels(['sampled training','baseline combat','same training','raw','different training'],fontsize=14)\n",
    "else:\n",
    "    ax.set_yticklabels(['baseline combat','same training','raw','sampled training','different training'],fontsize=14)\n",
    "ax.yaxis.set_label_position(\"right\")\n",
    "ax.yaxis.tick_right()\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/PearsonCorrelation_param{resolution}.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pearson.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pearson.std(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 0.6\n",
    "sns.set_style(\"whitegrid\")\n",
    "df_pearson = pd.read_csv(f'20210301_Pearson_correlation_param{resolution}.csv',index_col=0)\n",
    "ls_order = df_pearson.mean(axis=1).sort_values(ascending=False).index.tolist()\n",
    "%matplotlib inline\n",
    "s_date = '20201228'\n",
    "fig,ax = plt.subplots(figsize=(4.5,3),dpi=300)\n",
    "sns.boxplot(data=df_pearson.loc[ls_order].T,ax=ax,orient='h',showfliers=False,palette='muted')\n",
    "sns.stripplot(data=df_pearson.loc[ls_order].T,ax=ax,orient='h',palette='dark')\n",
    "ax.set_title(f'Cluster Correlation',fontsize=16) #\\n (resolution {resolution}\n",
    "ax.set_xlabel('Pearson Correlation',fontsize=14)\n",
    "if resolution == 0.6:\n",
    "    ax.set_yticklabels(['sampled training','baseline combat','same training','raw','different training'],fontsize=14)\n",
    "else:\n",
    "    pass\n",
    "    ax.set_yticklabels(['baseline combat','same training','sampled training','raw','different training'],fontsize=14)\n",
    "ax.yaxis.set_label_position(\"right\")\n",
    "ax.yaxis.tick_right()\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{codedir}/Figures/PearsonCorrelation_param.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pearson.mean(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pearson.std(axis=1)"
   ]
  },
  {
   "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": 4
}