Normalization_testing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalization of Adjacent TMA sections\n",
    "\n",
    "**Question:** What combination of transformation, scaling and batch regression produces the best batch correction on adjacent tissue sections stained with slighly different panel orders, but same antibody conjugates? \n",
    "\n",
    "**Samples:** \n",
    "- TMA: constructed from FFPE breast cancer cell lines and normal tonsil and breast tissue.\n",
    "- Sections: JE-TMA-41, JE-TMA-43, JE-TMA-62\n",
    "- Scenes: (i.e. TMA cores) \n",
    "  - 01, 07: Tonsil \n",
    "  - 04: Normal breast\n",
    "  - 02: HCC1143 (triple negative (TN), basal)\n",
    "  - 03: HCC3153 (TN, basal)\n",
    "  - 05, 06: T47D (luminal A, ER+, PR+)\n",
    "  - 08, 09: BT474 (luminal B, ER+, HER2+)\n",
    "  - 10, 11: AU565 (HER2+)\n",
    "  - 12, 13: MDAMB-436 (TN, claudin low)\n",
    "- Staining: The same fluorphore conjugated antibodies were applied, but not in the same order, and not the same lot.\n",
    "\n",
    "**Method**: We tested permutations of data transformation, scaling and batch regression. (see table) We evaluated effectiveness using the [kBET metric](https://github.com/theislab/kBET) to test for mixing of batches and the known subtype of the cell lines to ensure separation of biologically relevant cell types. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| Transformation | Scaling     | Regression  |\n",
    "|----------------|-------------|-------------|\n",
    "| raw            | standard    | regress_out |\n",
    "| log2           | min-max     | combat      |\n",
    "| arcsinh        | max-abs     |             |\n",
    "|                | robust      |             |\n",
    "|                | quantile    |             |\n",
    "|                | power       |             |\n",
    "|                | RESTORE     |             |\n",
    "|                | bg-quantile |             |"
   ]
  },
  {
   "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 importlib\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",
    "mpl.rc('figure', max_open_warning = 0)\n",
    "os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/')\n",
    "from mplex_image import visualize as viz, process, preprocess\n",
    "np.random.seed(1202)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#go to data location\n",
    "#change to correct directory\n",
    "os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cycIF_ValidationStudies/cycIF_Validation/Data')\n",
    "rootdir = os.getcwd()\n",
    "s_date = '20201228'\n",
    "if not os.path.exists(s_date):\n",
    "    os.mkdir(s_date)\n",
    "datadir = f'{rootdir}/validation_data'\n",
    "filterdir = f'{rootdir}/filtered_data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of contents <a name=\"contents\"></a>\n",
    "0. [Load Raw Data](#loadold)\n",
    "1. [Load Filtered Data](#load)\n",
    "2. [Scaling Function](#func)\n",
    "3. [Data Transformation](#run)\n",
    "4. [Scaling](#scale)\n",
    "   [Visualize](#norm)\n",
    "5. [kBET Results](#kbet)\n",
    "6. [Apply Batch Correction](#viznorm)\n",
    "8. [Cluster Composition](#cluster)\n",
    "9. [Combat fit/transform](#combat)\n",
    "9. [Combat fit/transform Application](#combatviz)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Raw Data  <a name=\"loadold\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(datadir)\n",
    "df_file = pd.DataFrame(index=os.listdir())\n",
    "df_file = df_file[df_file.index.str.contains('Filtered')]\n",
    "df_file['tissue'] = [item.split('_')[1] for item in df_file.index]\n",
    "df_file['dapi'] = ['DAPI' + item.split('y_DAPI')[1].split('.')[0] for item in df_file.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_sample = df_file.tissue.tolist()\n",
    "d_dapi = dict(zip(df_file.tissue.tolist(),df_file.dapi.tolist()))\n",
    "df_mi = pd.DataFrame()\n",
    "df_xy = pd.DataFrame()\n",
    "df_edge = pd.DataFrame()\n",
    "\n",
    "for s_sample in ls_sample:\n",
    "    print(f'loading {s_sample}')\n",
    "    df_mi = df_mi.append(pd.read_csv(f'{datadir}/features_{s_sample}_FilteredMeanIntensity_{d_dapi[s_sample]}.csv', index_col=0))\n",
    "    df_xy = df_xy.append(pd.read_csv(f'{datadir}/features_{s_sample}_CentroidXY.csv',index_col=0))\n",
    "    if os.path.exists(f'{datadir}/features_{s_sample}_EdgeCells153pixels_CentroidXY.csv'):\n",
    "        df_edge = df_edge.append(pd.read_csv(f'{datadir}/features_{s_sample}_EdgeCells153pixels_CentroidXY.csv',index_col=0))\n",
    "#scene naming\n",
    "df_mi.replace({'JE-TMA-43_scene13':'JE-TMA-43_scene12','JE-TMA-43_scene14':'JE-TMA-43_scene13'},inplace=True)\n",
    "df_mi['cell'] = [item.split('cell')[1] for item in df_mi.index]\n",
    "df_mi.index = df_mi.slide_scene + '_cell' + df_mi.cell \n",
    "\n",
    "df_xy.replace({'JE-TMA-43_scene13':'JE-TMA-43_scene12','JE-TMA-43_scene14':'JE-TMA-43_scene13'},inplace=True)\n",
    "df_xy['cell'] = [item.split('cell')[1] for item in df_xy.index]\n",
    "df_xy.index = df_xy.slide_scene + '_cell' + df_xy.cell "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select the right CD20 and CK5 fluorophore from JE-TMA-41\n",
    "ls_index = df_mi[df_mi.index.str.contains('JE-TMA-41') | df_mi.index.str.contains('SMT101Bx2-3')].index\n",
    "df_mi.loc[ls_index,'CD20_perinuc5'] = df_mi.loc[ls_index,'CD20P_perinuc5']\n",
    "df_mi.loc[ls_index,'CK5_cytoplasm'] = df_mi.loc[ls_index,'CK5P_cytoplasm']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_drop = ['AR_nuclei','HER2_cytoplasm','ER_nuclei25','CD44_nucadj2','Vim_perinuc5',\n",
    "           'DAPI1_nuclei','DAPI3_nuclei', 'DAPI4_nuclei',\n",
    "       'DAPI5_nuclei', 'DAPI6_nuclei', 'DAPI7_nuclei', 'DAPI8_nuclei',\n",
    "       'CSF1R_perinuc5', 'CAV1_perinuc5', 'BCL2_perinuc5', 'CoxIV_perinuc5',\n",
    "       'R0c2_perinuc5', 'R0c3_perinuc5', 'R0c4_perinuc5', 'R0c5_perinuc5',\n",
    "       'R5Qc2_perinuc5', 'R5Qc3_perinuc5', 'R5Qc4_perinuc5', 'R5Qc5_perinuc5',\n",
    "       'pS6RP_perinuc5', 'CK8_cytoplasm', 'EGFR_cytoplasm', 'MUC1_cytoplasm',\n",
    "       'DAPI0_nuclei', 'DAPI10_nuclei', 'DAPI5Q_nuclei', 'DAPI9_nuclei',\n",
    "       'FoxP3_nuclei', 'GRNZB_nuclei', 'H3K4_nuclei', 'PgR_nuclei',\n",
    "       'R0c2_nuclei', 'R0c3_nuclei', 'R0c4_nuclei', 'R0c5_nuclei',\n",
    "       'R5Qc2_nuclei', 'R5Qc3_nuclei', 'R5Qc4_nuclei', 'R5Qc5_nuclei',\n",
    "       'RAD51_nuclei', 'gH2AX_nuclei', 'pRB_nuclei', 'CD20P_perinuc5',\n",
    "       'CD8R_perinuc5', 'CD4R_perinuc5', 'ColI_perinuc5', 'ColIV_perinuc5',\n",
    "       'R8Qc2_perinuc5', 'R8Qc3_perinuc5', 'R8Qc4_perinuc5', 'R8Qc5_perinuc5',\n",
    "       'CK5P_cytoplasm', 'DAPI8Q_nuclei', 'LamB1_nuclei', 'R8Qc2_nuclei',\n",
    "       'R8Qc3_nuclei', 'R8Qc4_nuclei', 'R8Qc5_nuclei', 'BMP2_perinuc5',\n",
    "       'Glut1_perinuc5', 'PDGFRa_perinuc5', 'R1c2_perinuc5', 'pAKT_perinuc5',\n",
    "       'DAPI11_nuclei', 'DAPI12_nuclei', 'H3K27_nuclei', 'HIF1a_nuclei',\n",
    "       'LamB2_nuclei', 'R1c2_nuclei', 'cPARP_nuclei', 'pERK_nuclei']\n",
    "#drop markers not shared in all 3 TMAs\n",
    "df_filter_mi = df_mi.loc[:,~df_mi.columns.isin(ls_drop)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# filter out edge cores, edge cells in biopsies\n",
    "d_filter = {'JE-TMA-41_scene01':(df_xy.DAPI_Y > 5000),'JE-TMA-41_scene03':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-41_scene04':(df_xy.DAPI_Y < 1500),'JE-TMA-41_scene05':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-41_scene06':(df_xy.DAPI_Y < 1500),'JE-TMA-41_scene08':(df_xy.DAPI_Y < 1500),\n",
    "            'JE-TMA-41_scene09':(df_xy.DAPI_Y > 5000),'JE-TMA-41_scene11':(df_xy.DAPI_Y < 1500),\n",
    "            'JE-TMA-43_scene09':(df_xy.DAPI_Y < 1200),'JE-TMA-43_scene13':(df_xy.DAPI_Y < 1500),\n",
    "            #'JE-TMA-60_scene02':(df_xy.DAPI_Y > 5000),'JE-TMA-60_scene06':(df_xy.DAPI_Y < 1500),\n",
    "            #'JE-TMA-60_scene08':(df_xy.DAPI_Y > 5000),'JE-TMA-60_scene14':(df_xy.DAPI_X < 1500),\n",
    "            'JE-TMA-62_scene01':(df_xy.DAPI_Y > 5000),\n",
    "            'JE-TMA-62_scene02':(df_xy.DAPI_X > 5000),'JE-TMA-62_scene03':(df_xy.DAPI_X < 1000),\n",
    "            'JE-TMA-62_scene04':(df_xy.DAPI_Y < 1500),'JE-TMA-62_scene06':(df_xy.DAPI_X < 1000),\n",
    "            'JE-TMA-62_scene08':(df_xy.DAPI_Y > 5000),'JE-TMA-62_scene10':(df_xy.DAPI_Y < 1500),\n",
    "            'SMT101Bx2-3_scene002':(df_xy.DAPI_Y > 6000),'SMT101Bx2-5_scene002':(df_xy.DAPI_Y > 6000)}\n",
    "ls_filter_all = []\n",
    "for s_scene, filtercon in d_filter.items():\n",
    "    ls_filter = df_xy[(df_xy.slide_scene==s_scene) & filtercon].index.tolist()\n",
    "    ls_filter_all = ls_filter_all + ls_filter\n",
    "ls_filter = df_xy[(df_xy.slide_scene=='JE-TMA-60_scene02') &(df_xy.DAPI_X < 1500)].index.tolist()\n",
    "ls_filter_all = ls_filter_all + ls_filter\n",
    "#filter edge\n",
    "ls_filter_all = ls_filter_all + df_edge.index.tolist()\n",
    "df_filter_mi=df_filter_mi[~df_filter_mi.index.isin(ls_filter_all)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# check cells visually\n",
    "df_cluster = df_filter_mi.loc[:,['HER2_cellmem25','slide_scene']]\n",
    "df_cluster['cluster'] = 1\n",
    "df_cluster.drop('HER2_cellmem25',axis=1,inplace=True)\n",
    "\n",
    "viz.plot_clusters(df_cluster,df_xy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_filter_mi = df_filter_mi.loc[~df_filter_mi.index.isin(df_edge.index)]\n",
    "df_filter_mi.to_csv('20201228_SMT101-Bx2-Bx2-Bx4_FilteredMeanIntensity.csv')\n",
    "#df_filter_mi[df_filter_mi.index.str.contains('JE-TMA')].to_csv('20201210_JE-TMA-41-43-62_FilteredMeanIntensity.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Filtered Data <a name=\"load\"></a>\n",
    "\n",
    "Cells and markers filtered above\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(f'{filterdir}')\n",
    "df_mi = pd.read_csv('20201210_JE-TMA-41-43-62_FilteredMeanIntensity.csv',index_col=0)\n",
    "df_mi['batch'] = [item.split('_')[0] for item in df_mi.index]\n",
    "df_mi['scene'] = [item.split('_')[1] for item in df_mi.index]\n",
    "df_mi.drop('cell',axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Sampling <a name=\"run\"></a>\n",
    "\n",
    "Sample 2400 cells per TMA for kBET evaluation. Also calculate the bg quantile normalization.\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_sample = [#'SMT101Bx4-3','SMT101Bx2-3','SMT101Bx2-5',\n",
    " 'JE-TMA-41', 'JE-TMA-43', 'JE-TMA-62']\n",
    "ls_scene = ['scene01','scene07','scene04','scene03','scene05','scene08','scene11','scene12'] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(datadir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_bg = {}\n",
    "df_result = pd.DataFrame(index=[item.split('_')[0] for item in ls_marker],columns=ls_sample)\n",
    "for s_sample in ls_sample:\n",
    "    df_bg = pd.read_csv(f'features_{s_sample}_BackgroundQuantiles.csv',index_col=0)\n",
    "    se_bg = df_bg.loc[:,df_bg.columns.str.contains('_3')].mean()\n",
    "    se_bg.index = [item.split('_')[0] for item in se_bg.index]\n",
    "    for s_index in se_bg.index:\n",
    "        if df_result.index.isin([s_index]).sum()!=0:\n",
    "            df_result.loc[s_index,s_sample] = se_bg[s_index]\n",
    "    d_bg.update({s_sample:se_bg})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save\n",
    "if not os.path.exists('20201211_JE-TMA-41-43-62_bg-quant_normfactor.csv'):\n",
    "    df_result.dropna().to_csv('20201211_JE-TMA-41-43-62_bg-quant_normfactor.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(filterdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#select 1000 rand cells for each scene in ls_scene\n",
    "i_rand = 7\n",
    "df_sample = pd.DataFrame()\n",
    "for s_scene in ls_scene:\n",
    "    #print(s_scene)\n",
    "    df_scene = df_mi[df_mi.scene==s_scene]\n",
    "    for s_sample in ls_sample:\n",
    "        df_slide = df_scene[df_scene.batch==s_sample]\n",
    "        df_sample = df_sample.append(df_slide.sample(n=300,  replace=False,  random_state=i_rand)) #7, 8, 9 (2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save samples, bg quantile corrected\n",
    "df_bg_corrected = pd.DataFrame(index=df_sample.index)\n",
    "for s_sample in ls_sample:\n",
    "    df_bg = pd.read_csv(f'{datadir}/features_{s_sample}_BackgroundQuantiles.csv',index_col=0)\n",
    "    se_bg = df_bg.loc[:,df_bg.columns.str.contains('_3')].mean()\n",
    "    se_bg.index = [item.split('_')[0] for item in se_bg.index]\n",
    "    ls_index = df_sample[df_sample.batch==s_sample].index\n",
    "    for s_marker in df_sample.columns[df_sample.dtypes=='float64']:\n",
    "        if s_marker == 'DAPI2_nuclei':\n",
    "            continue\n",
    "        else:\n",
    "            df_bg_corrected.loc[ls_index,s_marker] =  df_sample.loc[ls_index,s_marker]/se_bg.loc[s_marker.split('_')[0]]\n",
    "df_bg_corrected['batch'] = [item.split('_')[0] for item in df_bg_corrected.index]\n",
    "df_bg_corrected.to_csv(f'2020120{i_rand}_JE-TMA-41-43-62_SampledMeanIntensity_raw_bg-quant.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot histograms\n",
    "%matplotlib inline\n",
    "s_trans = 'raw'\n",
    "bins=50\n",
    "for s_marker in df_bg_corrected.columns[df_bg_corrected.dtypes=='float64']:\n",
    "    print(s_marker)\n",
    "    fig,ax=plt.subplots(2,1,figsize = (5,5))\n",
    "    for idxs, s_batch in enumerate(sorted(set(df_sample.batch))):\n",
    "        df_batch = df_sample[(df_sample.batch==s_batch)].loc[:,s_marker] + 1 #set minimum to 1\n",
    "        if len(df_batch.dropna()) == 0:\n",
    "            continue\n",
    "        ax[0].hist(df_batch,bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "        ax[1].hist(df_bg_corrected.loc[df_bg_corrected.batch==s_batch,s_marker],bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "        ax[0].set_yscale('log')\n",
    "        ax[1].set_yscale('log')\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}: raw Unscaled Data')\n",
    "        ax[1].set_title(f'{s_marker.split(\"_\")[0]}: Subcellular BG')\n",
    "        ax[0].legend()\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{rootdir}/20201210/Different_Scaling_bg-quantile_{s_marker}_{s_trans}_{i_rand}.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save smpled data with transformations\n",
    "np.log2(df_sample.loc[:,df_sample.dtypes=='float64']+1).merge(df_sample.loc[:,['batch']],left_index=True,right_index=True).to_csv(f'2020120{i_rand}_JE-TMA-41-43-62_SampledMeanIntensity_log2.csv')\n",
    "np.arcsinh(df_sample.loc[:,df_sample.dtypes=='float64']).merge(df_sample.loc[:,['batch']],left_index=True,right_index=True).to_csv(f'2020120{i_rand}_JE-TMA-41-43-62_SampledMeanIntensity_arcsinh.csv')\n",
    "(df_sample.loc[:,df_sample.dtypes=='float64']).merge(df_sample.loc[:,['batch']],left_index=True,right_index=True).to_csv(f'2020120{i_rand}_JE-TMA-41-43-62_SampledMeanIntensity_raw.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling Function <a name=\"func\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from sklearn.preprocessing import minmax_scale\n",
    "from sklearn.preprocessing import MaxAbsScaler\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.preprocessing import RobustScaler\n",
    "from sklearn.preprocessing import Normalizer\n",
    "from sklearn.preprocessing import QuantileTransformer\n",
    "from sklearn.preprocessing import PowerTransformer\n",
    "#functions\n",
    "def test_scaling(X):\n",
    "    distributions = [\n",
    "     ('Unscaled data', X),\n",
    "     ('Data after standard scaling',\n",
    "        StandardScaler().fit_transform(X)),\n",
    "     ('Data after min-max scaling',\n",
    "        MinMaxScaler().fit_transform(X)),\n",
    "     ('Data after max-abs scaling',\n",
    "        MaxAbsScaler().fit_transform(X)),\n",
    "     ('Data after robust scaling',\n",
    "        RobustScaler(quantile_range=(25, 75),with_scaling=True).fit_transform(X)),\n",
    "     ('Data after power transformation (Box-Cox)',\n",
    "      PowerTransformer(method='box-cox').fit_transform(X)),\n",
    "     ('Data after quantile transformation (gaussian pdf)',\n",
    "        QuantileTransformer(output_distribution='normal')\n",
    "        .fit_transform(X)),\n",
    "      ]\n",
    "    return(distributions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Different Scaling <a name=\"scale\"></a>\n",
    "\n",
    "Test scaling methods: standard, robust, min-max, power, etc.\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#markers to norm\n",
    "ls_marker= ['HER2_cellmem25', 'Vim_nucadj2', 'CD20_perinuc5', 'CD3_perinuc5',\n",
    "       'CD31_perinuc5', 'CD4_perinuc5', 'CD44_perinuc5', 'CD45_perinuc5',\n",
    "       'CD68_perinuc5', 'CD8_perinuc5', 'PD1_perinuc5', 'PDPN_perinuc5',\n",
    "       'aSMA_perinuc5', 'CK14_cytoplasm', 'CK17_cytoplasm', 'CK19_cytoplasm',\n",
    "       'CK5_cytoplasm', 'CK7_cytoplasm', 'Ecad_cytoplasm', 'DAPI2_nuclei',\n",
    "       'ER_nuclei', 'Ki67_nuclei', 'LamAC_nuclei', 'PCNA_nuclei',\n",
    "       'pHH3_nuclei']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(ls_marker)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot all cores, different scaling\n",
    "%matplotlib inline\n",
    "s_trans = 'raw'#'arcsinh'#'log2'#\n",
    "bins=250\n",
    "for s_marker in ls_marker:\n",
    "    print(s_marker)\n",
    "    fig,ax=plt.subplots(5,1,figsize = (5,10))\n",
    "    for idxs, s_batch in enumerate(sorted(set(df_mi.batch))):\n",
    "        df_batch = df_mi[(df_mi.batch==s_batch)].loc[:,s_marker] + 1 #set minimum to 1\n",
    "        if len(df_batch.dropna()) == 0:\n",
    "            continue\n",
    "        if s_trans == 'log2':\n",
    "            X = np.log2(df_batch.values.reshape(-1,1)) + 1\n",
    "        elif s_trans == 'arcsinh':\n",
    "            X = np.arcsinh(df_batch.values.reshape(-1,1))\n",
    "        else:\n",
    "            X = df_batch.values.reshape(-1,1)\n",
    "        #print(X.min())\n",
    "        distributions = test_scaling(X)\n",
    "        for idx, tu_dist in enumerate(distributions):\n",
    "            ax[idx].hist(tu_dist[1],bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "            ax[idx].set_title(f'{s_trans} {tu_dist[0]}')\n",
    "            if s_trans == 'raw':\n",
    "                if idx == 5:\n",
    "                    pass\n",
    "                elif idx == 6:\n",
    "                    pass\n",
    "                else:\n",
    "                    ax[idx].set_yscale('log')\n",
    "            ax[idx].set_xlim(np.quantile(tu_dist[1],.001),np.quantile(tu_dist[1],.999))\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}: {s_trans} Unscaled Data')\n",
    "        ax[0].legend()\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{rootdir}/{s_date}/Different_Scaling_all_{s_marker}_{s_trans}.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load bg-quant and restore\n",
    "df_bgquant = pd.read_csv(f'{datadir}/20201211_JE-TMA-41-43-62_bg-quant_normfactor.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scaling compare standard, robust, restore and bg quantile\n",
    "%matplotlib inline\n",
    "s_trans = 'raw'#'arcsinh'#'log2'#\n",
    "bins=250\n",
    "for s_marker in ls_marker:\n",
    "    print(s_marker)\n",
    "    fig,ax=plt.subplots(7,1,figsize = (5,12))\n",
    "    for idxs, s_batch in enumerate(sorted(set(df_mi.batch))):\n",
    "        df_batch = df_mi[(df_mi.batch==s_batch)].loc[:,s_marker] + 1 #set minimum to 1\n",
    "        if len(df_batch.dropna()) == 0:\n",
    "            continue\n",
    "        if s_trans == 'log2':\n",
    "            X = np.log2(df_batch.values.reshape(-1,1)) + 1\n",
    "        else:\n",
    "            X = df_batch.values.reshape(-1,1)\n",
    "        dist_std = StandardScaler().fit_transform(X),\n",
    "        dist_robust = RobustScaler(quantile_range=(25, 75),with_scaling=True).fit_transform(X),\n",
    "        dist_restore = X/df_restore.loc[s_marker.split('_')[0],s_batch]\n",
    "        dist_bgquant = X/df_bgquant.loc[s_marker.split('_')[0],s_batch]\n",
    "        for idx, tu_dist in enumerate(distributions):\n",
    "            ax[idx].hist(tu_dist[1],bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "            ax[idx].set_title(f'{s_trans} {tu_dist[0]}')\n",
    "            if s_trans == 'raw':\n",
    "                ax[idx].set_yscale('log')\n",
    "            ax[idx].set_xlim(np.quantile(tu_dist[1],.001),np.quantile(tu_dist[1],.999))\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}: {s_trans} Unscaled Data')\n",
    "        ax[0].legend()\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{rootdir}/{s_date}/Different_Scaling_all_{s_marker}_{s_trans}.png')\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load sampled data,scale, save for kBET analysis (no bug)\n",
    "ls_sample = ['JE-TMA-41', 'JE-TMA-43', 'JE-TMA-62']\n",
    "i_rand = 9\n",
    "for s_trans in ['raw','log2','arcsinh']:\n",
    "    df_sample = pd.read_csv(f'{filterdir}/2020120{i_rand}_JE-TMA-41-43-62_SampledMeanIntensity_{s_trans}.csv',index_col=0)\n",
    "    #scaling results of sampled data\n",
    "    dd_result = {}\n",
    "    for s_batch in sorted(set(df_sample.batch)):\n",
    "        for s_marker in ls_marker:\n",
    "            d_result = {}\n",
    "            df_batch = df_sample.loc[df_sample.batch==s_batch,s_marker] + 1\n",
    "            X = df_batch.values.reshape(-1,1)\n",
    "            distributions = test_scaling(X)\n",
    "            for idx, tu_dist in enumerate(distributions[1:]):\n",
    "                d_result.update({tu_dist[0].split('after ')[1].split(' ')[0]:tu_dist[1]})\n",
    "            dd_result.update({f'{s_marker}_{s_batch}':d_result})\n",
    "    #convert to dataframe\n",
    "    df_result = pd.DataFrame()\n",
    "    for s_marker_batch, d_result in dd_result.items():\n",
    "        for s_norm, dist in d_result.items():\n",
    "            df_result[f'{s_marker_batch}_{s_norm}'] = dist.reshape(1,dist.shape[0])[0]\n",
    "    #save\n",
    "    ls_norm = sorted(set([item.split('_')[-1] for item in df_result.columns]))\n",
    "    for s_norm in ls_norm:\n",
    "        df_norm = df_result.loc[:,df_result.columns.str.contains(s_norm)]\n",
    "        df_out = pd.DataFrame()\n",
    "        for s_sample in ls_sample:\n",
    "            df_sample_norm = df_norm.loc[:,df_norm.columns.str.contains(s_sample)]\n",
    "            df_sample_norm.columns = [item.split('_')[0] for item in df_sample_norm.columns]\n",
    "            df_sample_norm['batch'] = s_sample\n",
    "            df_out = df_out.append(df_sample_norm)\n",
    "        df_out.index = df_sample.index\n",
    "        df_out.to_csv(f'{filterdir}/2020120{i_rand}_JE-TMA-41-43-62_SampledMeanIntensity_{s_trans}_{s_norm}.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#list files generated\n",
    "df_file = pd.DataFrame(index=os.listdir(filterdir))\n",
    "#df_file[df_file.index.str.contains('20201209_JE-TMA-41-43-62_SampledMeanIntensity')].index.tolist()\n",
    "df_file[df_file.index.str.contains('restore')].index.tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling Visualization  <a name=\"norm\"></a>\n",
    "\n",
    "The standard scaling method preformed the best in kbet analysis. Visualize histograms of results.\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#visualize\n",
    "\n",
    "#plot histograms\n",
    "#save\n",
    "df['batch'] = [item.split('_')[0] for item in df.index]\n",
    "#plot \n",
    "%matplotlib inline\n",
    "s_trans = ''\n",
    "s_date = \"all\"\n",
    "bins=50\n",
    "for s_marker in df.columns[df.dtypes=='float64']:\n",
    "    print(s_marker)\n",
    "    fig,ax=plt.subplots(2,1,figsize = (5,5))\n",
    "    for idxs, s_batch in enumerate(sorted(set(df.batch))):\n",
    "        df_batch = df[(df.batch==s_batch)].loc[:,s_marker] \n",
    "        if len(df_batch.dropna()) == 0:\n",
    "            continue\n",
    "        X = df_batch.values.reshape(-1,1)\n",
    "        dist_std = StandardScaler().fit_transform(X),\n",
    "        ax[0].hist(df.loc[df.index.str.contains(s_batch),s_marker],bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "        ax[1].hist(dist_std[0],bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "        ax[0].set_yscale('log')\n",
    "        ax[1].set_yscale('log')\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}: Raw Data')\n",
    "        ax[1].set_title(f'{s_marker.split(\"_\")[0]}: Standard Scaling')\n",
    "        ax[0].legend()\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{rootdir}/20201228/Different_Scaling_standard_{s_marker}_{s_trans}_{s_date}.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## kBET Results <a name=\"kbet\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "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('JE-TMA-41-43-62_kbet_')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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!='raw_restore']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_index= df_file.groupby('norm').mean_kbet.mean().sort_values().index\n",
    "df_file.groupby('norm').mean_kbet.mean().sort_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_file['Norm'] = pd.Categorical(\n",
    "    df_file['norm'], \n",
    "    categories=ls_index.tolist(), \n",
    "    ordered=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rootdir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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",
    "fig, ax = plt.subplots(figsize=(4,3.4),dpi=300)\n",
    "sns.lineplot(data=df_select.sort_values('mean_kbet'),x='norm',y='mean_kbet',ax=ax,err_style='bars')\n",
    "labels = [item.replace('div','global').replace('raw_','').replace('_','\\n').replace('raw','none') for item in df_select.groupby('norm').mean_kbet.mean().sort_values().index.tolist()]\n",
    "labels = ['z-score','combat','robust','regress\\nout','log2\\ncombat','RESTORE','raw']\n",
    "ax.set_xticks(range(len(df_select.groupby('norm').mean_kbet)))\n",
    "ax.set_xticklabels(labels,rotation=90)\n",
    "ax.set_ylabel('Rejection Rate')\n",
    "ax.set_xlabel('Normalization')\n",
    "ax.set_title('kBET Evaluation of\\n Batch Correction')\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "fig.savefig(f'{rootdir}/filtered_data/figures/BatchEffect_select.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "sns.set_style(\"whitegrid\")\n",
    "fig, ax = plt.subplots(figsize=(7,4),dpi=300)\n",
    "sns.lineplot(data=df_file.sort_values('Norm'),x='norm',y='mean_kbet',ax=ax,err_style='bars')\n",
    "labels = [item.replace('_out','').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=90)\n",
    "ax.set_ylabel('Rejection Rate')\n",
    "ax.set_xlabel('Normalization')\n",
    "ax.set_title('kBET Evaluation of Batch Correction I')\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "fig.savefig(f'{rootdir}/{s_date}/BatchEffectI.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"whitegrid\")\n",
    "#df_pearson = pd.read_csv(f'Pearson_correlation_all.csv',index_col=0)\n",
    "df_pearson = pd.read_csv(f'Pearson_correlation0.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,3),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",
    "ax.yaxis.set_label_position(\"right\")\n",
    "ax.set_yticklabels(['combat','log2 combat','RESTORE','regress out','raw','robust','z_score'],fontsize=14)\n",
    "ax.yaxis.tick_right()\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'{rootdir}/filtered_data/figures/PearsonCorrelation_all0.6.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "df_select = df_file[df_file.norm.isin(['raw_combat','raw_standard','restore_scale','raw','restore_div','restore_local','raw_regress_out'])] #'raw_restore_combat'\n",
    "sns.set_style(\"white\")\n",
    "fig, ax = plt.subplots(figsize=(4.2,3.2),dpi=300)\n",
    "sns.lineplot(data=df_select.sort_values('mean_kbet'),x='norm',y='mean_kbet',ax=ax,err_style='bars')\n",
    "labels = [item.replace('div','global').replace('raw_','').replace('_','\\n').replace('raw','none') for item in df_select.groupby('norm').mean_kbet.mean().sort_values().index.tolist()]\n",
    "ax.set_xticks(range(len(df_select.groupby('norm').mean_kbet)))\n",
    "ax.set_xticklabels(labels,rotation=90)\n",
    "ax.set_ylabel('Rejection Rate')\n",
    "ax.set_xlabel('Normalization')\n",
    "ax.set_title('kBET Evaluation of Batch Correction')\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "fig.savefig(f'{rootdir}/{s_date}/BatchEffect_select.png',dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_select.sort_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Batch Correction <a name=\"viznorm\"></a>\n",
    "\n",
    "Apply combat and regress out algorithims\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "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('20201207_JE-TMA-41-43-62_SampledMeanIntensity')]\n",
    "#df_file.index.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_index = [\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_robust.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_robust.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_combat.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_regress_out.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_power.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_min-max.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_max-abs.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_bg-quant.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_quantile.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2.csv',\n",
    " #'20201208_JE-TMA-41-43-62_SampledMeanIntensity_raw_restore.csv',\n",
    " #'20201209_JE-TMA-41-43-62_SampledMeanIntensity_raw_restore.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_quantile.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_max-abs.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_standard.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_power.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_min-max.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2_combat.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_regress_out.csv',\n",
    " #\n",
    " #top candidates\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_restore.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_combat.csv',\n",
    " '20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw_standard.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_restore_scale.csv',\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_restore_local.csv'\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_restore_div.csv'\n",
    " #'20201207_JE-TMA-41-43-62_SampledMeanIntensity_restore_div_arcsinh.csv'\n",
    " ]\n",
    "#ls_index = [ '20201208_JE-TMA-41-43-62_SampledMeanIntensity_raw.csv', '20201208_JE-TMA-41-43-62_SampledMeanIntensity_log2.csv',\n",
    "#            '20201207_JE-TMA-41-43-62_SampledMeanIntensity_raw.csv', '20201207_JE-TMA-41-43-62_SampledMeanIntensity_log2.csv',\n",
    "#            '20201209_JE-TMA-41-43-62_SampledMeanIntensity_raw.csv', '20201209_JE-TMA-41-43-62_SampledMeanIntensity_log2.csv']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_genes = [ 'CD45_perinuc5',  'CD20_perinuc5','PD1_perinuc5', 'CD3_perinuc5', 'CD4_perinuc5',\n",
    " 'CD8_perinuc5','CD68_perinuc5', 'aSMA_perinuc5','CD31_perinuc5', 'PDPN_perinuc5', 'Vim_nucadj2',\n",
    " 'CD44_perinuc5', 'CK14_cytoplasm', 'CK5_cytoplasm','CK17_cytoplasm',  'CK19_cytoplasm', 'CK7_cytoplasm', 'Ecad_cytoplasm',\n",
    " 'ER_nuclei','HER2_cellmem25', 'Ki67_nuclei', 'PCNA_nuclei', 'pHH3_nuclei', 'LamAC_nuclei']\n",
    "marker_genes = ['CD45', 'CD20', 'PD1', 'CD3', 'CD4', 'CD8', 'CD68', 'aSMA', 'CD31', 'PDPN', 'Vim', 'CD44',\n",
    "                'CK14','CK5', 'CK17', 'CK19', 'CK7', 'Ecad', 'ER', 'HER2', 'Ki67', 'PCNA', 'pHH3', 'LamAC']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "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",
    "for s_index in ls_index:\n",
    "    for s_norm in ['']: #,'none'#'regress_out','',,'combat'\n",
    "        print(s_index)\n",
    "        df = pd.read_csv(s_index.replace('kbet','SampledMeanIntensity'),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] for item in adata.obs.index]\n",
    "        adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]\n",
    "        #log transform, batch correct reduce dimensionality (PCA)\n",
    "        #if s_type=='Raw':\n",
    "        #    sc.pp.log1p(adata)\n",
    "        adata.raw = adata\n",
    "        #remove batch effect (non work well):#\n",
    "        print(s_norm)\n",
    "        if s_norm == 'combat':\n",
    "            sc.pp.combat(adata,key='batch')\n",
    "        elif s_norm == 'regress_out':\n",
    "            sc.pp.regress_out(adata, keys='batch')\n",
    "        elif s_norm == 'MNN':\n",
    "            sc.external.pp.mnn_correct(adata,batch_key='batch') #didn't work\n",
    "        else:\n",
    "            print('')\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=23)\n",
    "        sc.tl.umap(adata)\n",
    "        #umap plot\n",
    "        s_type = s_index.split('SampledMeanIntensity_')[1].split('.')[0]\n",
    "        figname = f\"UmapBatch_{s_type}_{s_norm}.png\"\n",
    "        fig,ax = plt.subplots(figsize=(4,4),dpi=300)\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=(4,4),dpi=300)\n",
    "        figname = f'UmapScene_{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",
    "        sc.tl.leiden(adata,resolution=0.25)\n",
    "        df = pd.DataFrame(data=adata.raw.X,index=adata.obs.index,columns=adata.var.index)\n",
    "        df['leiden'] = adata.obs['leiden']\n",
    "        #matrixplot\n",
    "        g=sns.clustermap(df.groupby('leiden').mean(),z_score=1,cmap='viridis')\n",
    "        ordered_genes=df.groupby('leiden').mean().iloc[:,g.dendrogram_col.reordered_ind].columns.tolist()\n",
    "        plt.close(g.fig)\n",
    "        sc.tl.dendrogram(adata,optimal_ordering =True, groupby=f'leiden')\n",
    "        fig,ax = plt.subplots(figsize=(8,4),dpi=300)\n",
    "        figname=f'Matrixplot_leiden_{s_type}_{s_norm}.png'\n",
    "        sc.pl.matrixplot(adata,  groupby=f'leiden',var_names=ordered_genes,\n",
    "                         dendrogram=True,ax=ax,save=figname,standard_scale='var',colorbar_title='Relative\\nintensity')\n",
    "        break\n",
    "    break\n",
    "sc.pl.pca_variance_ratio(adata, log=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cluster composition  <a name=\"cluster\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_je_tma = {'scene01':'tonsil1',\n",
    "    'scene02':'HCC1143', \n",
    "    'scene03':'HCC3153', \n",
    "    'scene04':'NBreast',\n",
    "    'scene05':'T47D',\n",
    "    'scene06':'T47D',\n",
    "    'scene07':'tonsil2',\n",
    "    'scene08':'BT474',\n",
    "    'scene09':'BT474',\n",
    "    'scene10':'AU565',\n",
    "    'scene11':'AU565',\n",
    "    'scene12':'MDAMB-436',\n",
    "    'scene13':'MDAMB-436',\n",
    " }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#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",
    "\n",
    "fig,ax=plt.subplots(figsize=(7,3.7), dpi=300)\n",
    "df_prop.columns = df_prop.columns.add_categories(['slide','scene'])\n",
    "df_prop.index = [item.replace('JE-TMA-','') 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,width=0.8)\n",
    "ax.legend(bbox_to_anchor=(1.02, 1.2), ncol=1)\n",
    "ax.set_ylabel('Fraction Positive')\n",
    "ax.set_title(f\"{s_trans.replace('_',' ')} {s_norm}\")\n",
    "ax.grid(False)\n",
    "plt.tight_layout()\n",
    "fig.savefig(f'./figures/StackedBar_{s_trans}_{s_norm}_Leiden.png')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# implement combat for fit/transform <a name=\"combat\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(f'{filterdir}')\n",
    "df = pd.read_csv('20201210_JE-TMA-41-43-62_FilteredMeanIntensity.csv',index_col=0)\n",
    "df['batch'] = [item.split('_')[0] for item in df.index]\n",
    "df['scene'] = [item.split('_')[1] for item in df.index]\n",
    "df.drop('cell',axis=1,inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load raw data\n",
    "#df = pd.read_csv(f'{s_index.replace(\".csv\",f\"_{s_norm}.csv\")}',index_col=0)\n",
    "#visualize raw\n",
    "adata = sc.AnnData(df.loc[df.index.str.contains('JE-TMA'),df.dtypes=='float64'])\n",
    "adata.obs['batch'] = [item.split('_scene')[0] for item in adata.obs.index]\n",
    "adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]\n",
    "sc.tl.pca(adata, svd_solver='auto')\n",
    "# calculate neighbors     \n",
    "sc.pp.neighbors(adata, n_neighbors=10, n_pcs=23)\n",
    "sc.tl.umap(adata)\n",
    "#umap plot\n",
    "sc.pl.umap(adata, color='batch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "#umap plot\n",
    "s_type = 'Full_TMA'\n",
    "s_norm = 'raw'\n",
    "figname = f\"UmapBatch_{s_type}_{s_norm}.png\"\n",
    "fig,ax = plt.subplots(figsize=(6,4))\n",
    "sc.pl.umap(adata, color='batch',title=f\"{s_type.replace('_',' ')} {s_norm}\",wspace=.25,ax=ax,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figname = f\"UmapScene_{s_type}_{s_norm}.png\"\n",
    "fig,ax = plt.subplots(figsize=(6,4))\n",
    "sc.pl.umap(adata, color='scene',title=f\"{s_type.replace('_',' ')} {s_norm}\",wspace=.25,ax=ax,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from https://github.com/brentp/combat.py/blob/master/combat.py\n",
    "import patsy\n",
    "import sys\n",
    "import numpy.linalg as la\n",
    "\n",
    "def aprior(gamma_hat):\n",
    "    m = gamma_hat.mean()\n",
    "    s2 = gamma_hat.var()\n",
    "    return (2 * s2 +m**2) / s2\n",
    "\n",
    "def bprior(gamma_hat):\n",
    "    m = gamma_hat.mean()\n",
    "    s2 = gamma_hat.var()\n",
    "    return (m*s2+m**3)/s2\n",
    "\n",
    "def it_sol(sdat, g_hat, d_hat, g_bar, t2, a, b, conv=0.0001):\n",
    "    n = (1 - np.isnan(sdat)).sum(axis=1)\n",
    "    g_old = g_hat.copy()\n",
    "    d_old = d_hat.copy()\n",
    "\n",
    "    change = 1\n",
    "    count = 0\n",
    "    while change > conv:\n",
    "        #print g_hat.shape, g_bar.shape, t2.shape\n",
    "        g_new = postmean(g_hat, g_bar, n, d_old, t2)\n",
    "        sum2 = ((sdat - np.dot(g_new.values.reshape((g_new.shape[0], 1)), np.ones((1, sdat.shape[1])))) ** 2).sum(axis=1)\n",
    "        d_new = postvar(sum2, n, a, b)\n",
    "       \n",
    "        change = max((abs(g_new - g_old) / g_old).max(), (abs(d_new - d_old) / d_old).max())\n",
    "        g_old = g_new #.copy()\n",
    "        d_old = d_new #.copy()\n",
    "        count = count + 1\n",
    "    adjust = (g_new, d_new)\n",
    "    return adjust \n",
    "\n",
    "def postmean(g_hat, g_bar, n, d_star, t2):\n",
    "    return (t2*n*g_hat+d_star * g_bar) / (t2*n+d_star)\n",
    "\n",
    "def postvar(sum2, n, a, b):\n",
    "    return (0.5 * sum2 + b) / (n / 2.0 + a - 1.0)\n",
    "\n",
    "def design_mat(mod, numerical_covariates, batch_levels):\n",
    "    # require levels to make sure they are in the same order as we use in the\n",
    "    # rest of the script.\n",
    "    design = patsy.dmatrix(\"~ 0 + C(batch, levels=%s)\" % str(batch_levels),\n",
    "                                                  mod, return_type=\"dataframe\")\n",
    "\n",
    "    mod = mod.drop([\"batch\"], axis=1)\n",
    "    numerical_covariates = list(numerical_covariates)\n",
    "    sys.stderr.write(\"found %i batches\\n\" % design.shape[1])\n",
    "    other_cols = [c for i, c in enumerate(mod.columns)\n",
    "                  if not i in numerical_covariates]\n",
    "    factor_matrix = mod[other_cols]\n",
    "    design = pd.concat((design, factor_matrix), axis=1)\n",
    "    if numerical_covariates is not None:\n",
    "        sys.stderr.write(\"found %i numerical covariates...\\n\"\n",
    "                            % len(numerical_covariates))\n",
    "        for i, nC in enumerate(numerical_covariates):\n",
    "            cname = mod.columns[nC]\n",
    "            sys.stderr.write(\"\\t{0}\\n\".format(cname))\n",
    "            design[cname] = mod[mod.columns[nC]]\n",
    "    sys.stderr.write(\"found %i categorical variables:\" % len(other_cols))\n",
    "    sys.stderr.write(\"\\t\" + \", \".join(other_cols) + '\\n')\n",
    "    return design\n",
    "\n",
    "def combat(data, batch, model=None, numerical_covariates=None):\n",
    "    \"\"\"Correct for batch effects in a dataset\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : pandas.DataFrame\n",
    "        A (n_features, n_samples) dataframe of the expression or methylation\n",
    "        data to batch correct\n",
    "    batch : pandas.Series\n",
    "        A column corresponding to the batches in the data, with index same as\n",
    "        the columns that appear in ``data``\n",
    "    model : patsy.design_info.DesignMatrix, optional\n",
    "        A model matrix describing metadata on the samples which could be\n",
    "        causing batch effects. If not provided, then will attempt to coarsely\n",
    "        correct just from the information provided in ``batch``\n",
    "    numerical_covariates : list-like\n",
    "        List of covariates in the model which are numerical, rather than\n",
    "        categorical\n",
    "    Returns\n",
    "    -------\n",
    "    corrected : pandas.DataFrame\n",
    "        A (n_features, n_samples) dataframe of the batch-corrected data\n",
    "    \"\"\"\n",
    "    if isinstance(numerical_covariates, str):\n",
    "        numerical_covariates = [numerical_covariates]\n",
    "    if numerical_covariates is None:\n",
    "        numerical_covariates = []\n",
    "\n",
    "    if model is not None and isinstance(model, pd.DataFrame):\n",
    "        model[\"batch\"] = list(batch)\n",
    "    else:\n",
    "        model = pd.DataFrame({'batch': batch})\n",
    "\n",
    "    batch_items = model.groupby(\"batch\").groups.items()\n",
    "    batch_levels = [k for k, v in batch_items]\n",
    "    batch_info = [v for k, v in batch_items]\n",
    "    n_batch = len(batch_info)\n",
    "    n_batches = np.array([len(v) for v in batch_info])\n",
    "    n_array = float(sum(n_batches))\n",
    "\n",
    "    # drop intercept\n",
    "    drop_cols = [cname for cname, inter in  ((model == 1).all()).iteritems() if inter == True]\n",
    "    drop_idxs = [list(model.columns).index(cdrop) for cdrop in drop_cols]\n",
    "    model = model[[c for c in model.columns if not c in drop_cols]]\n",
    "    numerical_covariates = [list(model.columns).index(c) if isinstance(c, str) else c\n",
    "            for c in numerical_covariates if not c in drop_cols]\n",
    "\n",
    "    design = design_mat(model, numerical_covariates, batch_levels)\n",
    "\n",
    "    sys.stderr.write(\"Standardizing Data across genes.\\n\")\n",
    "    #error shapes (3,7200) and (26,7200) not aligned: 7200 (dim 1) != 26 (dim 0)\n",
    "    B_hat = np.dot(np.dot(la.inv(np.dot(design.T, design)), design.T), data.T) #data.T\n",
    "    grand_mean = np.dot((n_batches / n_array).T, B_hat[:n_batch,:])\n",
    "    var_pooled = np.dot(((data - np.dot(design, B_hat).T)**2), np.ones((int(n_array), 1)) / int(n_array))\n",
    "\n",
    "    stand_mean = np.dot(grand_mean.T.reshape((len(grand_mean), 1)), np.ones((1, int(n_array))))\n",
    "    tmp = np.array(design.copy())\n",
    "    tmp[:,:n_batch] = 0\n",
    "    stand_mean  += np.dot(tmp, B_hat).T\n",
    "\n",
    "    s_data = ((data - stand_mean) / np.dot(np.sqrt(var_pooled), np.ones((1, int(n_array)))))\n",
    "\n",
    "    sys.stderr.write(\"Fitting L/S model and finding priors\\n\")\n",
    "    batch_design = design[design.columns[:n_batch]]\n",
    "    gamma_hat = np.dot(np.dot(la.inv(np.dot(batch_design.T, batch_design)), batch_design.T), s_data.T)\n",
    "\n",
    "    delta_hat = []\n",
    "\n",
    "    for i, batch_idxs in enumerate(batch_info):\n",
    "        #batches = [list(model.columns).index(b) for b in batches]\n",
    "        delta_hat.append(s_data[batch_idxs].var(axis=1))\n",
    "\n",
    "    gamma_bar = gamma_hat.mean(axis=1) \n",
    "    t2 = gamma_hat.var(axis=1)\n",
    "   \n",
    "\n",
    "    a_prior = list(map(aprior, delta_hat))\n",
    "    b_prior = list(map(bprior, delta_hat))\n",
    "\n",
    "    sys.stderr.write(\"Finding parametric adjustments\\n\")\n",
    "    gamma_star, delta_star = [], []\n",
    "    for i, batch_idxs in enumerate(batch_info):\n",
    "        #print '18 20 22 28 29 31 32 33 35 40 46'\n",
    "        #print batch_info[batch_id]\n",
    "\n",
    "        temp = it_sol(s_data[batch_idxs], gamma_hat[i],\n",
    "                     delta_hat[i], gamma_bar[i], t2[i], a_prior[i], b_prior[i])\n",
    "\n",
    "        gamma_star.append(temp[0])\n",
    "        delta_star.append(temp[1])\n",
    "\n",
    "    sys.stdout.write(\"Adjusting data\\n\")\n",
    "    bayesdata = s_data\n",
    "    gamma_star = np.array(gamma_star)\n",
    "    delta_star = np.array(delta_star)\n",
    "\n",
    "\n",
    "    for j, batch_idxs in enumerate(batch_info):\n",
    "\n",
    "        dsq = np.sqrt(delta_star[j,:])\n",
    "        dsq = dsq.reshape((len(dsq), 1))\n",
    "        denom =  np.dot(dsq, np.ones((1, n_batches[j])))\n",
    "        numer = np.array(bayesdata[batch_idxs] - np.dot(batch_design.loc[batch_idxs], gamma_star).T)\n",
    "\n",
    "        bayesdata[batch_idxs] = numer / denom\n",
    "   \n",
    "    vpsq = np.sqrt(var_pooled).reshape((len(var_pooled), 1))\n",
    "    bayesdata = bayesdata * np.dot(vpsq, np.ones((1, int(n_array)))) + stand_mean\n",
    " \n",
    "    return bayesdata\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#test function as is\n",
    "data = df.loc[:,df.dtypes=='float64'].T\n",
    "batch = df.batch\n",
    "bayesdata = combat(data, batch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#adapted from https://github.com/brentp/combat.py/blob/master/combat.py\n",
    "def combat_fit(data, batch, model=None, numerical_covariates=None):\n",
    "    \"\"\"Correct for batch effects in a dataset\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : pandas.DataFrame\n",
    "        A (n_features, n_samples) dataframe of the expression or methylation\n",
    "        data to batch correct\n",
    "    batch : pandas.Series\n",
    "        A column corresponding to the batches in the data, with index same as\n",
    "        the columns that appear in ``data``\n",
    "    model : patsy.design_info.DesignMatrix, optional\n",
    "        A model matrix describing metadata on the samples which could be\n",
    "        causing batch effects. If not provided, then will attempt to coarsely\n",
    "        correct just from the information provided in ``batch``\n",
    "    numerical_covariates : list-like\n",
    "        List of covariates in the model which are numerical, rather than\n",
    "        categorical\n",
    "    Returns\n",
    "    -------\n",
    "    gamma_star : centering parameters from combat fitting\n",
    "    delta_star : scaling parameters from combat fitting\n",
    "    \"\"\"\n",
    "    if isinstance(numerical_covariates, str):\n",
    "        numerical_covariates = [numerical_covariates]\n",
    "    if numerical_covariates is None:\n",
    "        numerical_covariates = []\n",
    "\n",
    "    if model is not None and isinstance(model, pd.DataFrame):\n",
    "        model[\"batch\"] = list(batch)\n",
    "    else:\n",
    "        model = pd.DataFrame({'batch': batch})\n",
    "\n",
    "    batch_items = model.groupby(\"batch\").groups.items()\n",
    "    batch_levels = [k for k, v in batch_items]\n",
    "    batch_info = [v for k, v in batch_items]\n",
    "    n_batch = len(batch_info)\n",
    "    n_batches = np.array([len(v) for v in batch_info])\n",
    "    n_array = float(sum(n_batches))\n",
    "\n",
    "    # drop intercept\n",
    "    drop_cols = [cname for cname, inter in  ((model == 1).all()).iteritems() if inter == True]\n",
    "    drop_idxs = [list(model.columns).index(cdrop) for cdrop in drop_cols]\n",
    "    model = model[[c for c in model.columns if not c in drop_cols]]\n",
    "    numerical_covariates = [list(model.columns).index(c) if isinstance(c, str) else c\n",
    "            for c in numerical_covariates if not c in drop_cols]\n",
    "\n",
    "    design = design_mat(model, numerical_covariates, batch_levels)\n",
    "\n",
    "    sys.stderr.write(\"Standardizing Data across genes.\\n\")\n",
    "    B_hat = np.dot(np.dot(la.inv(np.dot(design.T, design)), design.T), data.T) \n",
    "    grand_mean = np.dot((n_batches / n_array).T, B_hat[:n_batch,:])\n",
    "    var_pooled = np.dot(((data - np.dot(design, B_hat).T)**2), np.ones((int(n_array), 1)) / int(n_array))\n",
    "\n",
    "    stand_mean = np.dot(grand_mean.T.reshape((len(grand_mean), 1)), np.ones((1, int(n_array))))\n",
    "    tmp = np.array(design.copy())\n",
    "    tmp[:,:n_batch] = 0\n",
    "    stand_mean  += np.dot(tmp, B_hat).T\n",
    "\n",
    "    s_data = ((data - stand_mean) / np.dot(np.sqrt(var_pooled), np.ones((1, int(n_array)))))\n",
    "\n",
    "    sys.stderr.write(\"Fitting L/S model and finding priors\\n\")\n",
    "    batch_design = design[design.columns[:n_batch]]\n",
    "    gamma_hat = np.dot(np.dot(la.inv(np.dot(batch_design.T, batch_design)), batch_design.T), s_data.T)\n",
    "\n",
    "    delta_hat = []\n",
    "\n",
    "    for i, batch_idxs in enumerate(batch_info):\n",
    "        delta_hat.append(s_data[batch_idxs].var(axis=1))\n",
    "\n",
    "    gamma_bar = gamma_hat.mean(axis=1) \n",
    "    t2 = gamma_hat.var(axis=1)\n",
    "   \n",
    "\n",
    "    a_prior = list(map(aprior, delta_hat))\n",
    "    b_prior = list(map(bprior, delta_hat))\n",
    "\n",
    "    sys.stderr.write(\"Finding parametric adjustments\\n\")\n",
    "    gamma_star, delta_star = [], []\n",
    "    for i, batch_idxs in enumerate(batch_info):\n",
    "        temp = it_sol(s_data[batch_idxs], gamma_hat[i],\n",
    "                     delta_hat[i], gamma_bar[i], t2[i], a_prior[i], b_prior[i])\n",
    "\n",
    "        gamma_star.append(temp[0])\n",
    "        delta_star.append(temp[1])\n",
    "    return(gamma_star, delta_star)\n",
    "        \n",
    "def combat_transform(data, batch, gamma_star, delta_star,model=None, numerical_covariates=None):\n",
    "    \"\"\"Correct for batch effects in a dataset\n",
    "    Parameters\n",
    "    ----------\n",
    "    data : pandas.DataFrame\n",
    "        A (n_features, n_samples) dataframe of the expression or methylation\n",
    "        data to batch correct\n",
    "    batch : pandas.Series\n",
    "        A column corresponding to the batches in the data, with index same as\n",
    "        the columns that appear in ``data``\n",
    "    gamma_star : centering parameters from combat fitting\n",
    "    delta_star : scaling parameters from combat fitting\n",
    "    model : patsy.design_info.DesignMatrix, optional\n",
    "        A model matrix describing metadata on the samples which could be\n",
    "        causing batch effects. If not provided, then will attempt to coarsely\n",
    "        correct just from the information provided in ``batch``\n",
    "    numerical_covariates : list-like\n",
    "        List of covariates in the model which are numerical, rather than\n",
    "        categorical\n",
    "    Returns\n",
    "    -------\n",
    "    corrected : pandas.DataFrame\n",
    "        A (n_features, n_samples) dataframe of the batch-corrected data\n",
    "    \"\"\"\n",
    "    #get design\n",
    "    if isinstance(numerical_covariates, str):\n",
    "        numerical_covariates = [numerical_covariates]\n",
    "    if numerical_covariates is None:\n",
    "        numerical_covariates = []\n",
    "\n",
    "    if model is not None and isinstance(model, pd.DataFrame):\n",
    "        model[\"batch\"] = list(batch)\n",
    "    else:\n",
    "        model = pd.DataFrame({'batch': batch})\n",
    "    batch_items = model.groupby(\"batch\").groups.items()\n",
    "    batch_levels = [k for k, v in batch_items]\n",
    "    batch_info = [v for k, v in batch_items]\n",
    "    n_batch = len(batch_info)\n",
    "    n_batches = np.array([len(v) for v in batch_info])\n",
    "    n_array = float(sum(n_batches))\n",
    "    # drop intercept\n",
    "    drop_cols = [cname for cname, inter in  ((model == 1).all()).iteritems() if inter == True]\n",
    "    drop_idxs = [list(model.columns).index(cdrop) for cdrop in drop_cols]\n",
    "    model = model[[c for c in model.columns if not c in drop_cols]]\n",
    "    numerical_covariates = [list(model.columns).index(c) if isinstance(c, str) else c\n",
    "            for c in numerical_covariates if not c in drop_cols]\n",
    "    design = design_mat(model, numerical_covariates, batch_levels)\n",
    "    #standardize\n",
    "    sys.stderr.write(\"Standardizing Data across genes.\\n\")\n",
    "    B_hat = np.dot(np.dot(la.inv(np.dot(design.T, design)), design.T), data.T) \n",
    "    grand_mean = np.dot((n_batches / n_array).T, B_hat[:n_batch,:])\n",
    "    var_pooled = np.dot(((data - np.dot(design, B_hat).T)**2), np.ones((int(n_array), 1)) / int(n_array))\n",
    "\n",
    "    stand_mean = np.dot(grand_mean.T.reshape((len(grand_mean), 1)), np.ones((1, int(n_array))))\n",
    "    tmp = np.array(design.copy())\n",
    "    tmp[:,:n_batch] = 0\n",
    "    stand_mean  += np.dot(tmp, B_hat).T\n",
    "    s_data = ((data - stand_mean) / np.dot(np.sqrt(var_pooled), np.ones((1, int(n_array)))))\n",
    "    batch_design = design[design.columns[:n_batch]]\n",
    "    # adjust data\n",
    "    sys.stdout.write(\"Adjusting data\\n\")\n",
    "    bayesdata = s_data\n",
    "    gamma_star = np.array(gamma_star)\n",
    "    delta_star = np.array(delta_star)\n",
    "    #for each batch\n",
    "    for j, batch_idxs in enumerate(batch_info):\n",
    "\n",
    "        dsq = np.sqrt(delta_star[j,:])\n",
    "        dsq = dsq.reshape((len(dsq), 1))\n",
    "        denom =  np.dot(dsq, np.ones((1, n_batches[j]))) #divide by sqrt delta_star\n",
    "        numer = np.array(bayesdata[batch_idxs] - np.dot(batch_design.loc[batch_idxs], gamma_star).T) #subtract gamma_star\n",
    "\n",
    "        bayesdata[batch_idxs] = numer / denom\n",
    "    #multiply by square root of variance and add mean\n",
    "    vpsq = np.sqrt(var_pooled).reshape((len(var_pooled), 1))\n",
    "    bayesdata = bayesdata * np.dot(vpsq, np.ones((1, int(n_array)))) + stand_mean\n",
    "    return bayesdata\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fit\n",
    "gamma_star, delta_star = combat_fit(data, batch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma_star[0] #centering factor\n",
    "delta_star[0][0]  #scaling factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#transform\n",
    "bayesdata = combat_transform(data,batch,gamma_star, delta_star)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize combat histograms <a name=\"combatviz\"></a>\n",
    "\n",
    "[contents](#contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_norm = bayesdata.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#visualize\n",
    "\n",
    "#plot histograms\n",
    "#save\n",
    "df_norm['batch'] = [item.split('_')[0] for item in df_norm.index]\n",
    "#plot \n",
    "%matplotlib inline\n",
    "s_trans = ''\n",
    "s_date = \"all\"\n",
    "bins=50\n",
    "for s_marker in df_norm.columns[df_norm.dtypes=='float64']:\n",
    "    print(s_marker)\n",
    "    fig,ax=plt.subplots(2,1,figsize = (5,5))\n",
    "    for idxs, s_batch in enumerate(sorted(set(df_norm.batch))):\n",
    "        df_batch = df_norm[(df_norm.batch==s_batch)].loc[:,s_marker] #+ 1 #set minimum to 1\n",
    "        if len(df_batch.dropna()) == 0:\n",
    "            continue\n",
    "        ax[1].hist(df_batch,bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "        ax[0].hist(df.loc[df.index.str.contains(s_batch),s_marker],bins=bins,alpha=0.4, color=f'C{idxs}',label=s_batch)\n",
    "        ax[0].set_yscale('log')\n",
    "        ax[1].set_yscale('log')\n",
    "        ax[0].set_title(f'{s_marker.split(\"_\")[0]}: Raw Data')\n",
    "        ax[1].set_title(f'{s_marker.split(\"_\")[0]}: Combat')\n",
    "        ax[0].legend()\n",
    "    plt.tight_layout()\n",
    "    fig.savefig(f'{rootdir}/20201228/Different_Scaling_combat_{s_marker}_{s_trans}_{s_date}.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test on Biopsies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#test on biopsy samples\n",
    "df_test = pd.read_csv('20201228_SMT101-Bx2-Bx2-Bx4_FilteredMeanIntensity.csv',index_col=0)\n",
    "df_test.drop(['slide_scene','cell'], axis=1,inplace=True)\n",
    "df_test['batch'] = [item.split('_scene')[0] for item in df_test.index]\n",
    "df_sample = pd.DataFrame()\n",
    "df_control = pd.DataFrame() \n",
    "for s_batch in ['SMT101Bx2-3', 'SMT101Bx2-5', 'SMT101Bx4-3']:\n",
    "    df_sample = df_sample.append(df_test[df_test.batch==s_batch])\n",
    "for s_batch in ['JE-TMA-41', 'JE-TMA-43', 'JE-TMA-62']:\n",
    "    df_control = df_control.append(df_test[df_test.batch==s_batch])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fit\n",
    "df=df_control\n",
    "data = df.loc[:,df.dtypes=='float64'].T\n",
    "batch = df.batch\n",
    "gamma_star, delta_star = combat_fit(data, batch)\n",
    "#transform\n",
    "df=df_sample\n",
    "data = df.loc[:,df.dtypes=='float64'].T\n",
    "batch = df.batch\n",
    "bayesdata = combat_transform(data,batch,gamma_star, delta_star)\n",
    "df_norm = bayesdata.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#transform control\n",
    "\n",
    "df=df_control\n",
    "data = df.loc[:,df.dtypes=='float64'].T\n",
    "batch = df.batch\n",
    "bayesdata = combat_transform(data,batch,gamma_star, delta_star)\n",
    "df_norm_control = bayesdata.T\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#raw\n",
    "#load normalized data\n",
    "df = df_sample.append(df_control)\n",
    "#visualize\n",
    "adata = sc.AnnData(df.loc[:,df.dtypes=='float64'])\n",
    "adata.obs['batch'] = [item.split('_scene')[0] for item in adata.obs.index]\n",
    "adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]\n",
    "sc.tl.pca(adata, svd_solver='auto')\n",
    "# calculate neighbors     \n",
    "sc.pp.neighbors(adata, n_neighbors=10, n_pcs=23)\n",
    "sc.tl.umap(adata)\n",
    "#umap plot\n",
    "sc.pl.umap(adata, color='batch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#umap plot\n",
    "s_type = 'Biopsy & control'\n",
    "s_norm = 'raw'\n",
    "figname = f\"UmapBatch_{s_type}_{s_norm}.png\"\n",
    "fig,ax = plt.subplots(figsize=(6,4))\n",
    "sc.pl.umap(adata, color='batch',title=f\"{s_type.replace('_',' ')} {s_norm}\",wspace=.25,ax=ax,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load normalized data\n",
    "df = df_norm.append(df_norm_control)\n",
    "#visualize\n",
    "#adata = sc.AnnData(df.loc[:,ls_markers])\n",
    "adata = sc.AnnData(df.loc[:,df.dtypes=='float64'])\n",
    "adata.obs['batch'] = [item.split('_scene')[0] for item in adata.obs.index]\n",
    "adata.obs['scene'] = [item.split('_')[1] for item in adata.obs.index]\n",
    "sc.tl.pca(adata, svd_solver='auto')\n",
    "# calculate neighbors     \n",
    "sc.pp.neighbors(adata, n_neighbors=10, n_pcs=19)\n",
    "sc.tl.umap(adata)\n",
    "#umap plot\n",
    "sc.pl.umap(adata, color='batch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#umap plot\n",
    "s_type = 'Biopsy & control'#'Control'\n",
    "s_norm = 'combat'\n",
    "figname = f\"UmapBatch_{s_type}_{s_norm}.png\"\n",
    "fig,ax = plt.subplots(figsize=(6,4))\n",
    "sc.pl.umap(adata, color='batch',title=f\"{s_type.replace('_',' ')} {s_norm}\",wspace=.25,ax=ax,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cluster leiden\n",
    "sc.tl.leiden(adata, resolution = 0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "figname = f\"UmapLeiden_{s_type}_{s_norm}.png\"\n",
    "sc.pl.umap(adata, color='leiden',wspace=.25,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adata_smt = adata[adata.obs['leiden'].isin({'5', '6'}),:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(adata_smt, 'leiden', method='wilcoxon')\n",
    "sc.pl.rank_genes_groups(adata_smt, n_genes=10, sharey=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_genes = ['CD45_perinuc5','CD3_perinuc5','PDPN_perinuc5','CK5_cytoplasm',\n",
    "                'CK17_cytoplasm','Ki67_nuclei','CD44_perinuc5','CD8_perinuc5',\n",
    "               'CK19_cytoplasm','pHH3_nuclei','CD20_perinuc5','Vim_nucadj2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_markers = ['HER2_cellmem25', 'Vim_nucadj2',  'CD3_perinuc5',#\n",
    "       'CD31_perinuc5', 'CD4_perinuc5', 'CD44_perinuc5', #'\n",
    "       'CD68_perinuc5', 'CD8_perinuc5', 'PD1_perinuc5', 'PDPN_perinuc5',\n",
    "       'aSMA_perinuc5', 'CK14_cytoplasm', 'CK17_cytoplasm', #\n",
    "       'CK5_cytoplasm',  'Ecad_cytoplasm', 'DAPI2_nuclei',#,\n",
    "       'ER_nuclei', 'Ki67_nuclei', 'LamAC_nuclei', 'PCNA_nuclei',\n",
    "              'pHH3_nuclei','CD20_perinuc5','CD45_perinuc5','CK19_cytoplasm','CK7_cytoplasm']\n",
    "#umap plot\n",
    "s_type = 'Biopsy & control'#'Control'\n",
    "s_norm = 'combat'\n",
    "figname = f\"UmapMarkers_{s_type}_{s_norm}.png\"\n",
    "sc.pl.umap(adata, color=marker_genes,wspace=.25,save=figname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#CD20,CK19,CD45,CK7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(6,4))\n",
    "figname = f'UmapScene_{s_type}_{s_norm}.png'\n",
    "fig = sc.pl.umap(adata, color='scene',save=figname,ax=ax,title=f\"{s_type.replace('_',' ')} {s_norm}\")"
   ]
  },
  {
   "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
}