Normalization_testing_HER2-N75.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Batch Normalization on tissues\n",
    "\n",
    "\n",
    "**Question:** How do different batch normalization algorithms perform on adjacent section tissues stained with the same panel?\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 standard scaling, RESTORE, regress_out or combat normalization on three replicates of 1800 cells sampled from the data (600 cells per core). We calculated the kBET rejection rate, which is 1 for no mixing of batches, and 0 for perfect mixing. "
   ]
  },
  {
   "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\n",
    "from sklearn.metrics import silhouette_score\n",
    "import matplotlib as mpl\n",
    "mpl.rc('figure', max_open_warning = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cmIF')\n",
    "#from mplex_image import visualize as viz, process, preprocess\n",
    "np.random.seed(1211)"
   ]
  },
  {
   "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. [Functions](#func)\n",
    "3. [Data Transformation](#run)\n",
    "4. [Scaling](#scale)\n",
    "5. [kBET Results](#kbet)\n",
    "6. [Visualize Normalization](#normviz)\n",
    "7. [Apply Batch Correction](#viznorm)\n",
    "8. [Cluster](#cluster)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data\n",
    "#os.chdir('/home/groups/graylab_share/OMERO.rdsStore/engje/Data/cycIF_ValidationStudies/cycIF_Validation/Data')\n",
    "rootdir = os.getcwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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}/Data/filtered_data'\n",
    "#os.chdir(datadir)"
   ]
  },
  {
   "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('BM-Her2N75-15-17-18_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.isin(['raw_restore_combat','raw_restore_regress_out','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": [
    "%matplotlib inline\n",
    "fig, ax = plt.subplots(figsize=(7,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('_out','').replace('_','\\n') for item in ls_index.tolist()]\n",
    "ax.set_xticks(range(len(df_file.groupby('norm').mean_kbet)))\n",
    "ax.set_xticklabels(labels,rotation=45)\n",
    "ax.set_ylabel('Rejection Rate')\n",
    "ax.set_xlabel('Normalization')\n",
    "ax.set_title('kBET Evaluation of Batch Effect BM-HER2')\n",
    "fig.set_tight_layout(True)\n",
    "plt.tight_layout\n",
    "#fig.savefig(f'{rootdir}/{s_date}/BatchEffect.png',dpi=200)"
   ]
  },
  {
   "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": [
    "ls_index = [\n",
    " #'20201207_BM-Her2N75-15-17-18_SampledMeanIntensity_raw_restore.csv',\n",
    " #'20201207_BM-Her2N75-15-17-18_SampledMeanIntensity_raw.csv',\n",
    " '20201208_BM-Her2N75-15-17-18_SampledMeanIntensity_raw_restore.csv',\n",
    " '20201208_BM-Her2N75-15-17-18_SampledMeanIntensity_raw.csv',\n",
    " '20201209_BM-Her2N75-15-17-18_SampledMeanIntensity_raw_restore.csv',\n",
    " '20201209_BM-Her2N75-15-17-18_SampledMeanIntensity_raw.csv',\n",
    "]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalize \n",
    "%matplotlib inline\n",
    "sc.set_figure_params(scanpy=True, fontsize=14)\n",
    "for s_index in ls_index:\n",
    "    for s_norm in ['combat','regress_out','']: #,'none'\n",
    "        print(s_index)\n",
    "        df = pd.read_csv(s_index.replace('kbet','SampledMeanIntensity'),index_col=0)\n",
    "        print(df.index[0])\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",
    "        #sc.pl.pca(adata,color='batch')\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=(6,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=(6,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",
    "sc.pl.pca_variance_ratio(adata, log=True)"
   ]
  },
  {
   "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
}