JN4_Out of Distribution Prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Jupyter Notebook 4: Out of Distribution Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kang Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stderr",
     "text": "Using TensorFlow backend.\n"
    }
   ],
   "source": [
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "import sys\n",
    "import util_loss as ul\n",
    "#import the package to use\n",
    "import beta_vae\n",
    "import dhsic_vae\n",
    "from dentate_features import *\n",
    "from all_obs_linear_classifier_package import *\n",
    "import os,glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "['CD4 Naive T', 'Mk', 'B Activated', 'pDC', 'CD4 Memory T', 'T activated', 'B', 'CD14 Mono', 'DC', 'NK', 'CD16 Mono', 'CD8 T']\n"
    }
   ],
   "source": [
    "data = sc.read(\"./data/kang_seurat_normalized.h5ad\")\n",
    "data_train_full = sc.read(\"./data/kang_seurat_normalized_train.h5ad\")\n",
    "data_validate_full = sc.read(\"./data/kang_seurat_normalized_validate.h5ad\")\n",
    "cells = list(set(data_train_full.obs[\"cell_type\"]))\n",
    "print(cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cell in cells:\n",
    "\n",
    "    data_train_full_temp = data_train_full[-(data_train_full.obs[\"cell_type\"]==cell)]\n",
    "    data_validate_full_temp = data_validate_full[-(data_validate_full.obs[\"cell_type\"]==cell)]\n",
    "    print(data_train_full_temp.obs[\"cell_type\"].value_counts())\n",
    "\n",
    "    data_train_full_temp = ul.shuffle_adata(data_train_full_temp)\n",
    "    data_validate_full_temp = ul.shuffle_adata(data_validate_full_temp)\n",
    "    \n",
    "    #Declaring parameters\n",
    "    z = 5\n",
    "    al = 100\n",
    "    c = 500\n",
    "\n",
    "    mod_path1 = \"./models_kang_leave_one/latent\"+str(z)+\"_alpha\"+str(al)+\"_c\"+str(c)+\"_\"+cell\n",
    "    scg_model = beta_vae.VAEArithKeras(x_dimension= data_train_full.shape[1],z_dimension=z, \n",
    "                                          model_to_use =mod_path1,alpha=al,c_max=c)                  \n",
    "    scg_model.train(data_train_full_temp,validation_data=data_validate_full_temp,\n",
    "                    n_epochs=2,shuffle=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Manipulating Latent Space dimensions\n",
    "simulate_multiple_cell() recovers the dropped cell type by manipulating dimensions from every cell type.\n",
    "'''\n",
    "\n",
    "#Reloading pre-trained data\n",
    "from simulate_cell import *\n",
    "\n",
    "path = \"latent5_alpha20_c30_CD4 Naive T\"\n",
    "cell_to_drop = 'CD4 Naive T'\n",
    "scg_model = beta_vae.VAEArithKeras(x_dimension= data.shape[1],z_dimension=5,model_to_use=path,\n",
    "                                       alpha=20,c_max=30)\n",
    "scg_model.restore_model()\n",
    "\n",
    "data_temp = data[-(data.obs[\"cell_type\"]==cell_to_drop)]\n",
    "simulate_multiple_cell(path=path,data=data,model=scg_model,z_dim=5,feature=\"cell_type\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "<tf.Variable 'Variable_4:0' shape=() dtype=float32>\nModel: \"VAE\"\n_________________________________________________________________\nLayer (type)                 Output Shape              Param #   \n=================================================================\ninput (InputLayer)           (None, 14053)             0         \n_________________________________________________________________\nencoder (Model)              (None, 5)                 11896810  \n_________________________________________________________________\ndecoder (Model)              (None, 14053)             11906853  \n=================================================================\nTotal params: 23,803,663\nTrainable params: 23,797,263\nNon-trainable params: 6,400\n_________________________________________________________________\n<tf.Variable 'Variable_4:0' shape=() dtype=float32>\n"
    }
   ],
   "source": [
    "'''\n",
    "To limit the range of simulation\n",
    "It can range between mean value of the Source cell.\n",
    "Latent Representation plots can be used to limit the manipulation to a more local range to generate clearer plots.\n",
    "'''\n",
    "\n",
    "path = \"latent5_alpha20_c30_CD4 Naive T\"\n",
    "cell_to_drop = 'CD4 Naive T'\n",
    "scg_model = beta_vae.VAEArithKeras(x_dimension= data.shape[1],z_dimension=5,model_to_use=path,\n",
    "                                       alpha=20,c_max=30)\n",
    "scg_model.restore_model()\n",
    "\n",
    "data_temp = data[-(data.obs[\"cell_type\"]==cell_to_drop)]\n",
    "feature = \"cell_type\"\n",
    "cell = \"B\"\n",
    "z_dim = 5\n",
    "\n",
    "variable_names = data.var_names\n",
    "data_latent = scg_model.to_latent(data.X)\n",
    "latent_df = pd.DataFrame(data_latent)\n",
    "latent_df[feature] = list(data.obs[feature])\n",
    "latent_df[\"condition\"] = list(data.obs[\"condition\"])\n",
    "try:\n",
    "    os.makedirs(path+\"/gene_heatmaps/\")\n",
    "except OSError:\n",
    "    pass\n",
    "x_dim = data.shape[1]\n",
    "data_ast = latent_df[latent_df[feature]==cell]\n",
    "data_ast = data_ast[data_ast[\"condition\"]==\"CTRL\"]\n",
    "cell_one = data_ast.iloc[[0],[0,1,2,3,4]]\n",
    "\n",
    "for dim in range(z_dim):\n",
    "    '''\n",
    "    a and b are the ranges for simulation. Change it according to local requirements. \n",
    "    a = data_ast[dim].mean()\n",
    "    '''\n",
    "    a = min(data_latent[:,dim]) \n",
    "    b = max(data_latent[:,dim])\n",
    "    increment_range = np.arange(a,b,0.01)\n",
    "    result_array = np.empty((0, x_dim))\n",
    "    for inc in increment_range:\n",
    "            cell_latent = cell_one\n",
    "            #print(cell_latent)\n",
    "            #print(cell_latent.shape)\n",
    "            cell_latent.iloc[:,dim] = inc\n",
    "            cell_recon = scg_model.reconstruct(cell_latent)\n",
    "            result_array = np.append(result_array,cell_recon,axis=0)\n",
    "\n",
    "    result_adata = sc.AnnData(result_array, obs={\"inc_vals\":increment_range},var={\"var_names\":variable_names})\n",
    "    result_adata.write(path+\"/gene_heatmaps/\"+str(cell)+\"_\"+str(dim)+\".h5ad\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Creating Regression plots, UMAPS and PCA plots\n",
    "'''\n",
    "\n",
    "from seurat_umaps_reg_plots import *\n",
    "\n",
    "sc.tl.rank_genes_groups(data, groupby=\"cell_type\", method='t-test')\n",
    "x_dim = data.shape[1]\n",
    "\n",
    "path = \"latent5_alpha20_c30_CD4 Naive T\"\n",
    "cell_to_drop = 'CD4 Naive T'\n",
    "cells = list(set(data.obs[\"cell_type\"]))\n",
    "\n",
    "scg_model = beta_vae.VAEArithKeras(x_dimension= data.shape[1],z_dimension=5,model_to_use=path,\n",
    "                                       alpha=20,c_max=30)\n",
    "scg_model.restore_model()\n",
    "\n",
    "df_list = generate_simulated_reg_plots(path=path+\"/gene_heatmaps/\",\n",
    "                            actual_data=data,clust_typ = cell_to_drop,cells=cells)\n",
    "df = pd.DataFrame(df_list,columns=[\"name\",\"r_sq_all\",\"r_sq_100\"])\n",
    "df.to_csv(path+\"/reg_mean.csv\",index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Directly using PCA \n",
    "'''\n",
    "\n",
    "sim_data = sc.read(\"latent5_alpha20_c30_CD4 Naive T/gene_heatmaps/B_4.h5ad\")\n",
    "print(sim_data)\n",
    "generate_simulated_pca(path = \"latent5_alpha20_c30_CD4 Naive T/gene_heatmaps/\",actual_data = data,clust_typ = \"CD4 Naive T\",source_cell=\"B\",sim_data=sim_data,first_cell=\"B\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dentate Gyrus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sc.read(\"./data/dentate_gyrus_normalized.h5ad\")\n",
    "#data_train_full = sc.read(\"./data/dentate_gyrus_normalized_train.h5ad\")\n",
    "#data_validate_full = sc.read(\"./data/dentate_gyrus_normalized_validate.h5ad\")\n",
    "cells = list(set(data_train_full.obs[\"clusters\"]))\n",
    "print(cells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cell in cells:\n",
    "\n",
    "    data_train_full_temp = data_train_full[-(data_train_full.obs[\"clusters\"]==cell)]\n",
    "    data_validate_full_temp = data_validate_full[-(data_validate_full.obs[\"clusters\"]==cell)]\n",
    "    print(data_train_full_temp.obs[\"clusters\"].value_counts())\n",
    "\n",
    "    data_train_full_temp = ul.shuffle_adata(data_train_full_temp)\n",
    "    data_validate_full_temp = ul.shuffle_adata(data_validate_full_temp)\n",
    "    \n",
    "    #Declaring parameters\n",
    "    z = 5\n",
    "    al = 100\n",
    "    c = 30\n",
    "\n",
    "    mod_path1 = \"./models_dentate_leave_one/latent\"+str(z)+\"_alpha\"+str(al)+\"_c\"+str(c)+\"_\"+cell\n",
    "    scg_model = beta_vae.VAEArithKeras(x_dimension= data_train_full.shape[1],z_dimension=z, \n",
    "                                          model_to_use =mod_path1,alpha=al,c_max=c)                  \n",
    "    scg_model.train(data_train_full_temp,validation_data=data_validate_full_temp,\n",
    "                    n_epochs=2,shuffle=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Manipulating Latent Space dimensions\n",
    "'''\n",
    "\n",
    "#Reloading pre-trained data\n",
    "from simulate_cell import *\n",
    "\n",
    "path = \"latent5_alpha100_c30_Granule immature\"\n",
    "cell_to_drop = 'Granule immature'\n",
    "scg_model = beta_vae.VAEArithKeras(x_dimension= data.shape[1],z_dimension=5,model_to_use=path,\n",
    "                                       alpha=100,c_max=30)\n",
    "scg_model.restore_model()\n",
    "\n",
    "data_temp = data[-(data.obs[\"clusters\"]==cell_to_drop)]\n",
    "simulate_multiple_cell(path=path,data=data,model=scg_model,z_dim=5,feature=\"clusters\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "To limit the range of simulation\n",
    "It can range between mean value of the Source cell.\n",
    "Latent Representation plots can be used to limit the manipulation to a more local range to generate clearer plots.\n",
    "'''\n",
    "\n",
    "path = \"latent5_alpha100_c30_Granule immature\"\n",
    "cell_to_drop = 'Granule immature'\n",
    "scg_model = beta_vae.VAEArithKeras(x_dimension= data.shape[1],z_dimension=5,model_to_use=path,\n",
    "                                       alpha=100,c_max=500)\n",
    "scg_model.restore_model()\n",
    "\n",
    "data_temp = data[-(data.obs[\"clusters\"]==cell_to_drop)]\n",
    "feature = \"clusters\"\n",
    "cell = \"Neuroblast\"\n",
    "z_dim = 5\n",
    "\n",
    "variable_names = data.var_names\n",
    "data_latent = scg_model.to_latent(data.X)\n",
    "latent_df = pd.DataFrame(data_latent)\n",
    "latent_df[feature] = list(data.obs[feature])\n",
    "try:\n",
    "    os.makedirs(path+\"/gene_heatmaps/\")\n",
    "except OSError:\n",
    "    pass\n",
    "x_dim = data.shape[1]\n",
    "data_ast = latent_df[latent_df[feature]==cell]\n",
    "cell_one = data_ast.iloc[[0],[0,1,2,3,4]]\n",
    "\n",
    "for dim in range(z_dim):\n",
    "    '''\n",
    "    a and b are the ranges for simulation. Change it according to local requirements. \n",
    "    a = data_ast[dim].mean()\n",
    "    '''\n",
    "    a = min(data_latent[:,dim]) \n",
    "    b = max(data_latent[:,dim])\n",
    "    increment_range = np.arange(a,b,0.01)\n",
    "    result_array = np.empty((0, x_dim))\n",
    "    for inc in increment_range:\n",
    "            cell_latent = cell_one\n",
    "            #print(cell_latent)\n",
    "            #print(cell_latent.shape)\n",
    "            cell_latent.iloc[:,dim] = inc\n",
    "            cell_recon = scg_model.reconstruct(cell_latent)\n",
    "            result_array = np.append(result_array,cell_recon,axis=0)\n",
    "\n",
    "    result_adata = sc.AnnData(result_array, obs={\"inc_vals\":increment_range},var={\"var_names\":variable_names})\n",
    "    result_adata.write(path+\"/gene_heatmaps/\"+str(cell)+\"_\"+str(dim)+\".h5ad\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.tl.rank_genes_groups(data, groupby=\"clusters\", method='t-test')\n",
    "x_dim = data.shape[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Creating Regression plots, UMAPS and PCA plots\n",
    "'''\n",
    "\n",
    "from dentate_umaps_reg_plots import *\n",
    "\n",
    "path = \"latent5_alpha100_c30_Granule immature\"\n",
    "cell_to_drop = 'Granule immature'\n",
    "cells = list(set(data.obs[\"clusters\"]))\n",
    "\n",
    "scg_model = beta_vae.VAEArithKeras(x_dimension= data.shape[1],z_dimension=5,model_to_use=path,\n",
    "                                       alpha=50,c_max=30)\n",
    "scg_model.restore_model()\n",
    "\n",
    "df_list = generate_simulated_reg_plots(path=path+\"/gene_heatmaps/\",\n",
    "                            actual_data=data,clust_typ = cell_to_drop,cells=cells)\n",
    "df = pd.DataFrame(df_list,columns=[\"name\",\"r_sq_all\",\"r_sq_100\"])\n",
    "df.to_csv(path+\"/gene_heatmaps/reg_mean.csv\",index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.5-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}