genie3_evaluation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nethome/mbafna6/miniconda3/envs/CeSpGRN/lib/python3.8/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import networkx as nx\n",
    "import plotly.io as pio\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "import igviz as ig\n",
    "from node2vec import Node2Vec\n",
    "from gensim.models import KeyedVectors\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "\n",
    "# import src.preprocess as pre\n",
    "# import src.visualize as vis\n",
    "# pio.renderers.default = \"png\"\n",
    "\n",
    "from src import models, training\n",
    "from torch_geometric.nn import GAE\n",
    "import torch\n",
    "import os\n",
    "import genie3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_cells = 1597"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"./data/scMultiSim_final/scmultisim_dataframe.csv\")\n",
    "df5 = df[df['Cell_Type'] == \"2_1_5\"].iloc[:,4:].values\n",
    "df4 = df[df['Cell_Type'] == \"2_1_4\"].iloc[:,4:].values\n",
    "df3 = df[df['Cell_Type'] == \"2_1_3\"].iloc[:,4:].values\n",
    "df2 = df[df['Cell_Type'] == \"2_1_2\"].iloc[:,4:].values\n",
    "df1 = df[df['Cell_Type'] == \"2_1_1\"].iloc[:,4:].values\n",
    "\n",
    "countlist = [df1,df2,df3,df4,df5]\n",
    "gene_names = df.columns[4:]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cell Type 1 GRN, counts shape (240, 120)\n",
      "Tree method: RF, K: sqrt, Number of trees: 1000 \n",
      "\n",
      "Cell Type 2 GRN, counts shape (240, 120)\n",
      "Tree method: RF, K: sqrt, Number of trees: 1000 \n",
      "\n",
      "Cell Type 3 GRN, counts shape (240, 120)\n",
      "Tree method: RF, K: sqrt, Number of trees: 1000 \n",
      "\n",
      "Cell Type 4 GRN, counts shape (240, 120)\n",
      "Tree method: RF, K: sqrt, Number of trees: 1000 \n",
      "\n",
      "Cell Type 5 GRN, counts shape (240, 120)\n",
      "Tree method: RF, K: sqrt, Number of trees: 1000 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "celltype_grns = []\n",
    "for i,counts in enumerate(countlist):\n",
    "    print(f\"Cell Type {i+1} GRN, counts shape {counts.shape}\")\n",
    "    type_grn = genie3.GENIE3(counts, gene_names=gene_names, nthreads=8)\n",
    "    np.save(f\"./benchmark/genie3/scmultisim_final/CellType{i+1}GRN.npy\", type_grn)\n",
    "    celltype_grns.append(type_grn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_grn_path = \"./data/scMultiSim_final/simulated/cellgrn/\"\n",
    "\n",
    "gt_grns = np.zeros((len(df),len(gene_names),60))\n",
    "\n",
    "for i, filename in enumerate(os.listdir(gt_grn_path)):\n",
    "    grn = pd.read_csv(os.path.join(gt_grn_path, filename)).values\n",
    "    \n",
    "    gt_grns[i, : ,:] = grn\n",
    "    \n",
    "gt_grns.shape, gt_grns\n",
    "np.save(\"./data/scMultiSim_final/gt_grns.npy\", gt_grns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1200, 120, 60)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gt_grns = np.load(\"./data/scMultiSim_final/gt_grns.npy\")\n",
    "gt_grns.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "genie3_grns = np.zeros((len(df),len(gene_names),60))\n",
    "\n",
    "for i, row in df.iterrows():\n",
    "    type_num = row[\"Cell_Type\"].split(\"_\")[-1]\n",
    "    celltype_grn = np.load(f\"./benchmark/genie3/scmultisim_final/CellType{type_num}GRN.npy\")[:,:60]\n",
    "    genie3_grns[i,:,:] = celltype_grn\n",
    "np.save(\"./benchmark/genie3/scmultisim_final/genie3_grns.npy\", genie3_grns)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "gt_grns = np.load(\"./data/scMultiSim_final/gt_grns.npy\")\n",
    "genie3_grns = np.load(\"./benchmark/genie3/scmultisim_final/genie3_grns.npy\")\n",
    "cespgrn_grns = np.load(\"./out/scmultisim_final/1_preprocessing_output/initial_grns.npy\")\n",
    "clarify_grns = np.load(\"./benchmark/genie3/clarify_grns.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.01853888888888889, 1.4029852482084098)"
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.metrics import roc_auc_score, jaccard_score, accuracy_score, f1_score, precision_recall_curve, auc\n",
    "\n",
    "total_genes = 1200 * 120 * 60\n",
    "\n",
    "network_density = np.sum(gt_grns.flatten()>0) / total_genes\n",
    "\n",
    "p, r, _ = precision_recall_curve(gt_grns.flatten()>0, genie3_grns.flatten())\n",
    "network_density, auc(r,p) / network_density\n",
    "\n",
    "# auc(r,p) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.01853888888888889, 1.474749717495629)"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# total_genes = 1200 * 60 * 59\n",
    "# network_density = np.sum(gt_grns[:,:45,:45].flatten()) / total_genes\n",
    "\n",
    "\n",
    "p, r, _ = precision_recall_curve(gt_grns[:,:45,:45].flatten()>0, clarify_grns.flatten())\n",
    "network_density, auc(r,p) / network_density\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.026814814814814812"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "f1_score(gt_grns[:,:45,:45].flatten()>4,clarify_grns.flatten()>0, average=\"micro\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "inverted_gt = (gt_grns[:,:45,:45].flatten()>4) == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.026814814814814812"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f1_score(inverted_gt,(clarify_grns.flatten()>0)==0, average=\"micro\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "module 'numpy' has no attribute 'rand'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn [110], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m roc_auc_score(gt_grns\u001b[39m.\u001b[39mflatten()\u001b[39m>\u001b[39m\u001b[39m0\u001b[39m, np\u001b[39m.\u001b[39;49mrand\u001b[39m.\u001b[39mpermute(gt_grns\u001b[39m.\u001b[39mflatten()\u001b[39m>\u001b[39m\u001b[39m0\u001b[39m))\n",
      "File \u001b[0;32m~/miniconda3/envs/CeSpGRN/lib/python3.8/site-packages/numpy/__init__.py:311\u001b[0m, in \u001b[0;36m__getattr__\u001b[0;34m(attr)\u001b[0m\n\u001b[1;32m    308\u001b[0m     \u001b[39mfrom\u001b[39;00m \u001b[39m.\u001b[39;00m\u001b[39mtesting\u001b[39;00m \u001b[39mimport\u001b[39;00m Tester\n\u001b[1;32m    309\u001b[0m     \u001b[39mreturn\u001b[39;00m Tester\n\u001b[0;32m--> 311\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mAttributeError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mmodule \u001b[39m\u001b[39m{!r}\u001b[39;00m\u001b[39m has no attribute \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m    312\u001b[0m                      \u001b[39m\"\u001b[39m\u001b[39m{!r}\u001b[39;00m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39mformat(\u001b[39m__name__\u001b[39m, attr))\n",
      "\u001b[0;31mAttributeError\u001b[0m: module 'numpy' has no attribute 'rand'"
     ]
    }
   ],
   "source": [
    "roc_auc_score(gt_grns.flatten()>0, np.rand.permute(gt_grns.flatten()>0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "CeSpGRN",
   "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.8.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "a3478565159e01bf18e9a76bba0cc0584a9ca455954cded0bbc47b7868cb1082"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}