Vignette.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chronos Vignette"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This vignette walks through a simple exercise in training Chronos on a subset of DepMap public 20Q4 and the Sanger Institute's Project Score data. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import chronos\n",
    "import os\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "from taigapy import default_tc as tc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some tweaks that will make plots more legible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import rcParams\n",
    "rcParams['axes.titlesize'] = 14\n",
    "rcParams['axes.spines.right'] = False\n",
    "rcParams['axes.spines.top'] = False\n",
    "rcParams['savefig.dpi'] = 200\n",
    "rcParams['savefig.transparent'] = False\n",
    "rcParams['font.family'] = 'Arial'\n",
    "rcParams['font.size'] = '11'\n",
    "rcParams['figure.dpi'] = 200\n",
    "rcParams[\"savefig.facecolor\"] = (1, 1, 1.0, 0.2)\n",
    "\n",
    "rcParams['xtick.labelsize'] = 10\n",
    "rcParams['ytick.labelsize'] = 10\n",
    "rcParams['legend.fontsize'] = 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chronos always requires at least three dataframes: \n",
    "* a matrix of readcounts with sequenced entities as the index, individual sgRNAs as the columns, and values indicating how many reads were found for that sgRNA. A sequenced entity any vector of sgRNA readcounts read out during the experiment. It could be a sequencing run of pDNA, or of a biological replicate at some time point during the experiment.\n",
    "* A sequence map mapping sequenced entities to either pDNA or a cell line and giving the days since infection and pDNA batch. \n",
    "* A guide map mapping sgRNAs to genes. Each sgRNA included must map to one and only one gene."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we'll load a small subset of the DepMap Avana data. The files have been reformatted from the release to the format Chronos expects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sequence_map = pd.read_csv(\"Data/SampleData/AvanaSequenceMap.csv\")\n",
    "guide_map = pd.read_csv(\"Data/SampleData/AvanaGuideMap.csv\")\n",
    "readcounts = chronos.read_hdf5(\"Data/SampleData/AvanaReadcounts.hdf5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sequence maps must have the columns\n",
    "\n",
    "* sequence_id (str), which must match a row in readcounts\n",
    "* cell_line_name (str). Must be \"pDNA\" for pDNA, and each pDNA batch must have at least one pDNA measurement.\n",
    "* pDNA batch (any simple hashable type, preferably int or str). pDNA measurements sharing the same batch will be grouped and averaged, then used as the reference for all biological replicate sequencings assigned that same batch. If you don't have multiple pDNA batches (by far the most common experimental condition), just fill this column with 0 or some other constant value.\n",
    "* days: days post infection. This value will be ignored for pDNA.\n",
    "\n",
    "Other columns will be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sequence_ID</th>\n",
       "      <th>ScreenID</th>\n",
       "      <th>days</th>\n",
       "      <th>pDNA_batch</th>\n",
       "      <th>Replicate</th>\n",
       "      <th>ScreenType</th>\n",
       "      <th>cell_line_name</th>\n",
       "      <th>ModelConditionID</th>\n",
       "      <th>Library</th>\n",
       "      <th>PassesQC</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>HEL-311Cas9_RepA_p4_Avana-3</td>\n",
       "      <td>SC-000004.AV01</td>\n",
       "      <td>21</td>\n",
       "      <td>Avana-3</td>\n",
       "      <td>A</td>\n",
       "      <td>2DS</td>\n",
       "      <td>ACH-000004</td>\n",
       "      <td>MC-000004-pA3k</td>\n",
       "      <td>Avana</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>HEL-311Cas9_RepB_p4_Avana-3</td>\n",
       "      <td>SC-000004.AV01</td>\n",
       "      <td>21</td>\n",
       "      <td>Avana-3</td>\n",
       "      <td>B</td>\n",
       "      <td>2DS</td>\n",
       "      <td>ACH-000004</td>\n",
       "      <td>MC-000004-pA3k</td>\n",
       "      <td>Avana</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>KU812-311cas9-RepA-p6_Avana-3</td>\n",
       "      <td>SC-000074.AV01</td>\n",
       "      <td>21</td>\n",
       "      <td>Avana-3</td>\n",
       "      <td>A</td>\n",
       "      <td>2DS</td>\n",
       "      <td>ACH-000074</td>\n",
       "      <td>MC-000074-OKtM</td>\n",
       "      <td>Avana</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>KU812-311cas9-RepB-p6_Avana-3</td>\n",
       "      <td>SC-000074.AV01</td>\n",
       "      <td>21</td>\n",
       "      <td>Avana-3</td>\n",
       "      <td>B</td>\n",
       "      <td>2DS</td>\n",
       "      <td>ACH-000074</td>\n",
       "      <td>MC-000074-OKtM</td>\n",
       "      <td>Avana</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>T47D-311Cas9-RepA-p6_Avana-4</td>\n",
       "      <td>SC-000147.AV01</td>\n",
       "      <td>21</td>\n",
       "      <td>Avana-4</td>\n",
       "      <td>A</td>\n",
       "      <td>2DS</td>\n",
       "      <td>ACH-000147</td>\n",
       "      <td>MC-000147-Uovr</td>\n",
       "      <td>Avana</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     sequence_ID        ScreenID  days pDNA_batch Replicate  \\\n",
       "0    HEL-311Cas9_RepA_p4_Avana-3  SC-000004.AV01    21    Avana-3         A   \n",
       "1    HEL-311Cas9_RepB_p4_Avana-3  SC-000004.AV01    21    Avana-3         B   \n",
       "2  KU812-311cas9-RepA-p6_Avana-3  SC-000074.AV01    21    Avana-3         A   \n",
       "3  KU812-311cas9-RepB-p6_Avana-3  SC-000074.AV01    21    Avana-3         B   \n",
       "4   T47D-311Cas9-RepA-p6_Avana-4  SC-000147.AV01    21    Avana-4         A   \n",
       "\n",
       "  ScreenType cell_line_name ModelConditionID Library  PassesQC  \n",
       "0        2DS     ACH-000004   MC-000004-pA3k   Avana      True  \n",
       "1        2DS     ACH-000004   MC-000004-pA3k   Avana      True  \n",
       "2        2DS     ACH-000074   MC-000074-OKtM   Avana      True  \n",
       "3        2DS     ACH-000074   MC-000074-OKtM   Avana      True  \n",
       "4        2DS     ACH-000147   MC-000147-Uovr   Avana      True  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sequence_map[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Guide maps must have the columns \n",
    "\n",
    "* sgrna (str): must match a column in readcounts. An sgrna can only appear once in this data frame.\n",
    "* gene (str): the gene the sgrna maps to.\n",
    "\n",
    "Other columns will be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sgrna</th>\n",
       "      <th>GenomeAlignment</th>\n",
       "      <th>gene</th>\n",
       "      <th>nAlignments</th>\n",
       "      <th>DropReason</th>\n",
       "      <th>UsedByChronos</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>AAAAATGCGCAAATTCAGCG</td>\n",
       "      <td>chr3_138742712_-</td>\n",
       "      <td>PIK3CB (5291)</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>AAAACACATCAGTATAACAT</td>\n",
       "      <td>chr3_49368469_+</td>\n",
       "      <td>RHOA (387)</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>AAAACTACAGAAGCCTCCCG</td>\n",
       "      <td>chr10_34450424_-</td>\n",
       "      <td>PARD3 (56288)</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>AAAAGGCCTGACATATCTGA</td>\n",
       "      <td>chr15_66444677_+</td>\n",
       "      <td>MAP2K1 (5604)</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  sgrna   GenomeAlignment           gene  nAlignments  \\\n",
       "0  AAAAATGCGCAAATTCAGCG  chr3_138742712_-  PIK3CB (5291)          1.0   \n",
       "1  AAAACACATCAGTATAACAT   chr3_49368469_+     RHOA (387)          1.0   \n",
       "2  AAAACTACAGAAGCCTCCCG  chr10_34450424_-  PARD3 (56288)          1.0   \n",
       "3  AAAAGGCCTGACATATCTGA  chr15_66444677_+  MAP2K1 (5604)          2.0   \n",
       "\n",
       "   DropReason  UsedByChronos  \n",
       "0         NaN           True  \n",
       "1         NaN           True  \n",
       "2         NaN           True  \n",
       "3         NaN           True  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "guide_map[:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, here's what readcounts should look like. They can include NaNs. Note the axes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>AAAAATGCGCAAATTCAGCG</th>\n",
       "      <th>AAAACACATCAGTATAACAT</th>\n",
       "      <th>AAAACTACAGAAGCCTCCCG</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>HEL-311Cas9_RepA_p4_Avana-3</th>\n",
       "      <td>101.0</td>\n",
       "      <td>224.0</td>\n",
       "      <td>636.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>HEL-311Cas9_RepB_p4_Avana-3</th>\n",
       "      <td>147.0</td>\n",
       "      <td>400.0</td>\n",
       "      <td>350.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>KU812-311cas9-RepA-p6_Avana-3</th>\n",
       "      <td>124.0</td>\n",
       "      <td>191.0</td>\n",
       "      <td>364.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>KU812-311cas9-RepB-p6_Avana-3</th>\n",
       "      <td>129.0</td>\n",
       "      <td>536.0</td>\n",
       "      <td>1280.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                               AAAAATGCGCAAATTCAGCG  AAAACACATCAGTATAACAT  \\\n",
       "HEL-311Cas9_RepA_p4_Avana-3                   101.0                 224.0   \n",
       "HEL-311Cas9_RepB_p4_Avana-3                   147.0                 400.0   \n",
       "KU812-311cas9-RepA-p6_Avana-3                 124.0                 191.0   \n",
       "KU812-311cas9-RepB-p6_Avana-3                 129.0                 536.0   \n",
       "\n",
       "                               AAAACTACAGAAGCCTCCCG  \n",
       "HEL-311Cas9_RepA_p4_Avana-3                   636.0  \n",
       "HEL-311Cas9_RepB_p4_Avana-3                   350.0  \n",
       "KU812-311cas9-RepA-p6_Avana-3                 364.0  \n",
       "KU812-311cas9-RepB-p6_Avana-3                1280.0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "readcounts.iloc[:4, :3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To QC the data, we'll want control groups. We'll use predefined sets of common and nonessential genes, and use these to define control sets of sgRNAs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "common_essentials = pd.read_csv(\"Data/SampleData/AchillesCommonEssentialControls.csv\")[\"Gene\"]\n",
    "nonessentials = pd.read_csv(\"Data/SampleData/AchillesNonessentialControls.csv\")[\"Gene\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "positive_controls = guide_map.sgrna[guide_map.gene.isin(common_essentials)]\n",
    "negative_controls = guide_map.sgrna[guide_map.gene.isin(nonessentials)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalizing Readcounts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's recommended that you normalize the readcounts before running `nan_outgrowths` below to correct for library size/sequencing depth effects. Chronos' `normalize_readcounts` function has different behavior depending on what is provided. In the ideal case that both a `sequence_map` and abundant negative controls (>=1000) are supplied, normalization process in two stages:\n",
    "- pDNA counts are scaled so that the modes of the readcount log distributions are aligned.\n",
    "- Each other sequence is then scaled such that the median of the log reads of negative controls within the sequence matches the median of the same guides in the negative controls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "readcounts = chronos.normalize_readcounts(readcounts, negative_controls, sequence_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NaNing clonal outgrowths"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Achilles, we've observed rare instances where a single guide in a single biological replicate will produce an unexpectedly large number of readcounts, while other guides targeting the same gene or other replicates of the same cell line do not show many readcounts. We suspect this is the result of a single clone gaining some fitness advantage. Although it _could_ be related to a change induced by the guide, in general it's probably misleading. Therefore Chronos has an option to identify and remove these events. Note this change occurs in place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating LFC\n",
      "stacking and annotating LFC\n",
      "finding group boundaries\n",
      "removing cases with only one guide and replicate\n",
      "finding maximal values\n",
      "making mask\n",
      "found 219 outgrowths, 6.9E-04 of the total\n",
      "pivoting mask\n",
      "aligning_mask\n",
      "NaNing\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/model.py:476: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  lfc_stack[\"Mask\"].loc[bad_rows] = True\n"
     ]
    }
   ],
   "source": [
    "chronos.nan_outgrowths(readcounts=readcounts, guide_gene_map=guide_map,\n",
    "                                   sequence_map=sequence_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QCing the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can generate a report with basic QC metrics about your data. You don't have to have control guides to do this, but the report is most useful if you do. If you don't have the `reportlab` python package installed, this section will error and should be skipped. This command will write a pdf report named \"Initial QC.pdf\" in the `./Data/reports` directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "reportdir = \"./Data/reports\"\n",
    "# permanently deletes the directory - careful if you edit this line!\n",
    "! rm -rf \"./Data/reports\"\n",
    "! mkdir \"./Data/reports\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/homebrew/lib/python3.11/site-packages/umap/distances.py:1063: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  @numba.jit()\n",
      "/opt/homebrew/lib/python3.11/site-packages/umap/distances.py:1071: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  @numba.jit()\n",
      "/opt/homebrew/lib/python3.11/site-packages/umap/distances.py:1086: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  @numba.jit()\n",
      "/opt/homebrew/lib/python3.11/site-packages/umap/umap_.py:660: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.\n",
      "  @numba.jit()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculating replicate correlation\n",
      "generating control separation metrics\n",
      "Plotting log fold-change distribution\n",
      "plotting control separation metrics\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/evaluations.py:317: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n",
      "/Users/dempster/packages/chronos_public/chronos/reports.py:194: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n"
     ]
    }
   ],
   "source": [
    "from chronos import reports\n",
    "metrics = reports.qc_initial_data(\"Initial QC\", readcounts, sequence_map,guide_map, \n",
    "        negative_controls, positive_controls,\n",
    "                  directory=reportdir\n",
    "       )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look in the Data/reports directory to see the QC report, \"Initial QC.pdf\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train Chronos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we initialize the model. Note the form of the data: each of the three parameters is actually a dictionary. If we were training the model with data from multiple libraries simultaneously, each library's data would have its own entries in the dict. \n",
    "\n",
    "The `negative_control_sgrnas` is an optional parameter, but including it will allow 1. better removal of library size effects from readcounts, and 2. estimation of the negative binomial quadratic overdispersion parameter per screen, which is otherwise a fixed hyperparameter. If provided, these should be cutting sgRNAs that are strongly expected to have no viability impact.\n",
    "\n",
    "`log_dir` is an optional argument containing a directory for tensorflow to write summaries to. We include it here so that tensorboard can load the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "logdir = \"./Data/logs\"\n",
    "# permanently deletes the directory - careful if you edit this line!\n",
    "! rm -rf \"./Data/logs\"\n",
    "! mkdir \"./Data/logs\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normalizing readcounts\n",
      "\n",
      "\n",
      "Finding all unique guides and genes\n",
      "found 3474 unique guides and 883 unique genes in avana\n",
      "found 3474 unique guides and 883 unique genes overall\n",
      "\n",
      "finding guide-gene mapping indices\n",
      "\n",
      "finding all unique sequenced replicates, cell lines, and pDNA batches\n",
      "found 92 unique sequences (excluding pDNA) and 44 unique cell lines in avana\n",
      "found 92 unique replicates and 44 unique cell lines overall\n",
      "\n",
      "finding replicate-cell line mappings indices\n",
      "\n",
      "finding replicate-pDNA mappings indices\n",
      "\n",
      "\n",
      "assigning float constants\n",
      "Estimating or aligning variances\n",
      "\tEstimating excess variance (alpha) for avana\n",
      "Creating excess variance tensors\n",
      "\tCreated excess variance tensor for avana with shape [92, 1]\n",
      "initializing graph\n",
      "\n",
      "building gene effect mask\n",
      "\n",
      "building doubling vectors\n",
      "made days vector of shape [92, 1] for avana\n",
      "\n",
      "building late observed timepoints\n",
      "\tbuilt normalized timepoints for avana with shape (92, 3474) (replicates X guides)\n",
      "\n",
      "building t0 reads\n",
      "\n",
      "\n",
      "Building variables\n",
      "\n",
      "building t0 reads estimate\n",
      "made t0 batch with shape (3, 3474) for avana\n",
      "building guide efficacy\n",
      "built guide efficacy: shape [1, 3475]\n",
      "building growth rate\n",
      "built growth rate: shape {'avana': [44, 1]}\n",
      "\n",
      "building cell line efficacy\n",
      "built cell line efficacy: shapes {'avana': [44, 1]}\n",
      "building screen delay\n",
      "built screen delay\n",
      "building gene effect\n",
      "built core gene effect: 44 cell lines by 883 genes\n",
      "\n",
      "\n",
      "Connecting graph nodes in model\n",
      "\n",
      "building effective days\n",
      "built effective days, shapes {'avana': [92, 883]}\n",
      "\n",
      "building gene effect growth graph nodes\n",
      "built gene effect growth graph nodes, shapes {'avana': [44, 883]}\n",
      "\n",
      "building combined efficacy\n",
      "built combined efficacy, shape {'avana': [44, 3474]}\n",
      "built expanded combined efficacy, shapes {'avana': [92, 3474]}\n",
      "\n",
      "building growth estimates of edited cells and overall estimates of fold change in guide abundance\n",
      "built growth and change\n",
      "\n",
      "building unnormalized estimates of final abundance\n",
      "built unnormalized abundance\n",
      "\n",
      "building normalized estimates of final abundance\n",
      "built normalized abundance\n",
      "initializing precost variables\n",
      "\n",
      "\n",
      "Building all costs\n",
      "\n",
      "assembling guide efficacy regularization\n",
      "building smoothed regularization\n",
      "\n",
      "building t0 reads regularization/cost\n",
      "\n",
      "building NB2 cost\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/model.py:377: UserWarning: Between 0 (batch=Index(['avana_Avana-2', 'avana_Avana-3', 'avana_Avana-4'], dtype='object')) and 0 (batch=Index(['avana_Avana-2', 'avana_Avana-3', 'avana_Avana-4'], dtype='object')) negative control sgRNAs were found to be systematically over- or under-represented in the screens and excluded.\n",
      "  warn(\"Between %i (batch=%r) and %i (batch=%r) negative control sgRNAs were found to be \\\n",
      "2024-06-26 13:56:19.005308: W tensorflow/c/c_api.cc:304] Operation '{name:'excess_variance/avana/Assign' id:7 op device:{requested: '', assigned: ''} def:{{{node excess_variance/avana/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](excess_variance/avana, excess_variance/avana/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building other regularizations\n",
      "\n",
      "Creating optimizer\n",
      "\tcreating log at ./Data/logs\n",
      "initializing rest of graph\n",
      "estimating initial screen efficacy and gene effect\n",
      "\t avana\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:56:19.429765: W tensorflow/c/c_api.cc:304] Operation '{name:'GE/library_effect/avana/Assign' id:150 op device:{requested: '', assigned: ''} def:{{{node GE/library_effect/avana/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](GE/library_effect/avana, GE/library_effect/avana/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\tverifying graph integrity\n",
      "verifying user inputs\n",
      "verifying variables\n",
      "verifying calculated terms\n",
      "\tavana _gene_effect\n",
      "\tavana _selected_efficacies\n",
      "\tavana_predicted_readcounts_unscaled\n",
      "\tavana _predicted_readcounts\n",
      "\tavana _normalized_readcounts\n",
      "\tavana _cost_presum\n",
      "sess run\n",
      "\tavana _cost\n",
      "\tavana _full_costs\n",
      "ready to train\n"
     ]
    }
   ],
   "source": [
    "model = chronos.Chronos(\n",
    "    sequence_map={\"avana\": sequence_map},\n",
    "    guide_gene_map={\"avana\": guide_map},\n",
    "    readcounts={\"avana\": readcounts},\n",
    "    negative_control_sgrnas={\"avana\": negative_controls},\n",
    "    log_dir=logdir\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have tensorboard, the cell below will show Chronos' node structure. `GE` means gene effect (relative change in growth rate), `FC` means predicted fold change, `t0` is the inferred relative guide abundance at t0, and `out_norm` is the predicted readcounts. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "      <iframe id=\"tensorboard-frame-afb3aa6ac55279a8\" width=\"100%\" height=\"800\" frameborder=\"0\">\n",
       "      </iframe>\n",
       "      <script>\n",
       "        (function() {\n",
       "          const frame = document.getElementById(\"tensorboard-frame-afb3aa6ac55279a8\");\n",
       "          const url = new URL(\"/\", window.location);\n",
       "          const port = 6006;\n",
       "          if (port) {\n",
       "            url.port = port;\n",
       "          }\n",
       "          frame.src = url;\n",
       "        })();\n",
       "      </script>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%reload_ext tensorboard\n",
    "!kill $(ps -e | grep 'tensorboard' | awk '{print $1}')\n",
    "%tensorboard --logdir ./data/logs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, optimizing the model:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we train  the model for 301 epochs. This should take a minute or so with periodic updates provided"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NB2 cost 0.3381509608150151\n",
      "Full cost 0.4240511122616703\n",
      "relative_growth_rate\n",
      "\tavana max 1.043, min 0.95978\n",
      "mean guide efficacy 0.9992352771862419\n",
      "t0_offset SD: [('avana', 7.063544027959644e-05)]\n",
      "\n",
      "gene mean -0.21254807454822527\n",
      "SD of gene means 0.2991141928308989\n",
      "Mean of gene SDs 0.2578709824610634\n",
      "\n",
      "\n",
      "\n",
      "51 epochs trained, time taken 0:00:02, projected remaining 0:00:08\n",
      "NB2 cost 0.2240884994186334\n",
      "Full cost 0.2644521088205495\n",
      "relative_growth_rate\n",
      "\tavana max 1.414, min 0.55608\n",
      "mean guide efficacy 0.9392223699161444\n",
      "t0_offset SD: [('avana', 0.1440585015866918)]\n",
      "\n",
      "gene mean -0.036717414513631544\n",
      "SD of gene means 0.43679826354833046\n",
      "Mean of gene SDs 0.2226218661144294\n",
      "\n",
      "\n",
      "\n",
      "101 epochs trained, time taken 0:00:03, projected remaining 0:00:06\n",
      "NB2 cost 0.21280469953893333\n",
      "Full cost 0.23817233417840725\n",
      "relative_growth_rate\n",
      "\tavana max 1.588, min 0.53601\n",
      "mean guide efficacy 0.8823429599704864\n",
      "t0_offset SD: [('avana', 0.13702414530460086)]\n",
      "\n",
      "gene mean 0.0002514219071480029\n",
      "SD of gene means 0.44750184303968366\n",
      "Mean of gene SDs 0.19451889618482882\n",
      "\n",
      "\n",
      "\n",
      "151 epochs trained, time taken 0:00:05, projected remaining 0:00:04\n",
      "NB2 cost 0.21201140651746014\n",
      "Full cost 0.23698614499778542\n",
      "relative_growth_rate\n",
      "\tavana max 1.556, min 0.52784\n",
      "mean guide efficacy 0.8777429933323182\n",
      "t0_offset SD: [('avana', 0.13656130436660804)]\n",
      "\n",
      "gene mean -0.0008748366306366722\n",
      "SD of gene means 0.45079206466631006\n",
      "Mean of gene SDs 0.1894554426614766\n",
      "\n",
      "\n",
      "\n",
      "201 epochs trained, time taken 0:00:06, projected remaining 0:00:03\n",
      "NB2 cost 0.2124745700852519\n",
      "Full cost 0.2371299701941426\n",
      "relative_growth_rate\n",
      "\tavana max 1.517, min 0.52762\n",
      "mean guide efficacy 0.8716317776295865\n",
      "t0_offset SD: [('avana', 0.1353208222639242)]\n",
      "\n",
      "gene mean 0.00018050873232808103\n",
      "SD of gene means 0.4526064138884017\n",
      "Mean of gene SDs 0.18727941704491946\n",
      "\n",
      "\n",
      "\n",
      "251 epochs trained, time taken 0:00:08, projected remaining 0:00:01\n",
      "NB2 cost 0.21192968418406305\n",
      "Full cost 0.23662779407818568\n",
      "relative_growth_rate\n",
      "\tavana max 1.531, min 0.52565\n",
      "mean guide efficacy 0.8697248417595845\n",
      "t0_offset SD: [('avana', 0.13544495870550516)]\n",
      "\n",
      "gene mean 6.740398555972341e-05\n",
      "SD of gene means 0.4551613563147576\n",
      "Mean of gene SDs 0.18580156362472813\n",
      "\n",
      "\n",
      "\n",
      "301 epochs trained, time taken 0:00:09, projected remaining 0:00:00\n",
      "NB2 cost 0.21219287011369709\n",
      "Full cost 0.23663274152108432\n",
      "relative_growth_rate\n",
      "\tavana max 1.519, min 0.53036\n",
      "mean guide efficacy 0.8658095869613612\n",
      "t0_offset SD: [('avana', 0.13438702495566443)]\n",
      "\n",
      "gene mean -0.0002363982701545945\n",
      "SD of gene means 0.455523036612303\n",
      "Mean of gene SDs 0.18558964501833553\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model.train(301, report_freq=50, burn_in_period=50, ge_only=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## After Training"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving and Restoring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chronos' `save` method dumps all the inputs, outputs, and model parameters to the specified directory. These files are written such that they can be read in individually and analyzed, but also used to restore the model by passing the directory path to the function `load_saved_model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "savedir = \"Data/Achilles_run\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isdir(savedir):\n",
    "    os.mkdir(savedir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.save(savedir, overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved files:\n",
      "\n",
      "\tlibrary_effect.csv\n",
      "\tcell_line_growth_rate.csv\n",
      "\tavana_predicted_readcounts.hdf5\n",
      "\tparameters.json\n",
      "\tguide_efficacy.csv\n",
      "\tavana_sequence_map.csv\n",
      "\tt0_offset.csv\n",
      "\tavana_predicted_lfc.hdf5\n",
      "\tavana_guide_gene_map.csv\n",
      "\tgene_effect_corrected.hdf5\n",
      "\tavana_negative_control_sgrnas.csv\n",
      "\tscreen_delay.csv\n",
      "\tscreen_excess_variance.csv\n",
      "\tcell_line_efficacy.csv\n",
      "\tgene_effect.hdf5\n",
      "\tavana_readcounts.hdf5\n"
     ]
    }
   ],
   "source": [
    "print(\"Saved files:\\n\\n\" + '\\n'.join(['\\t' + s for s in os.listdir(savedir)\n",
    "                if s.endswith(\"csv\")\n",
    "                or s.endswith(\"hdf5\")\n",
    "                or s.endswith(\"json\")\n",
    "                ]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The .hdf5 files are binaries written with chronos' `write_hdf5` function, which is an efficient method for writing large matrices. They can be read with chronos' `read_hdf5` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Restoring the model can be done with a single function call:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Finding all unique guides and genes\n",
      "found 3474 unique guides and 883 unique genes in avana\n",
      "found 3474 unique guides and 883 unique genes overall\n",
      "\n",
      "finding guide-gene mapping indices\n",
      "\n",
      "finding all unique sequenced replicates, cell lines, and pDNA batches\n",
      "found 92 unique sequences (excluding pDNA) and 44 unique cell lines in avana\n",
      "found 92 unique replicates and 44 unique cell lines overall\n",
      "\n",
      "finding replicate-cell line mappings indices\n",
      "\n",
      "finding replicate-pDNA mappings indices\n",
      "\n",
      "\n",
      "assigning float constants\n",
      "Estimating or aligning variances\n",
      "\tEstimating excess variance (alpha) for avana\n",
      "Creating excess variance tensors\n",
      "\tCreated excess variance tensor for avana with shape [92, 1]\n",
      "initializing graph\n",
      "\n",
      "building gene effect mask\n",
      "\n",
      "building doubling vectors\n",
      "made days vector of shape [92, 1] for avana\n",
      "\n",
      "building late observed timepoints\n",
      "\tbuilt normalized timepoints for avana with shape (92, 3474) (replicates X guides)\n",
      "\n",
      "building t0 reads\n",
      "\n",
      "\n",
      "Building variables\n",
      "\n",
      "building t0 reads estimate\n",
      "made t0 batch with shape (3, 3474) for avana\n",
      "building guide efficacy\n",
      "built guide efficacy: shape [1, 3475]\n",
      "building growth rate\n",
      "built growth rate: shape {'avana': [44, 1]}\n",
      "\n",
      "building cell line efficacy\n",
      "built cell line efficacy: shapes {'avana': [44, 1]}\n",
      "building screen delay\n",
      "built screen delay\n",
      "building gene effect\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/model.py:377: UserWarning: Between 0 (batch=Index(['avana_avana_Avana-2', 'avana_avana_Avana-3', 'avana_avana_Avana-4'], dtype='object')) and 0 (batch=Index(['avana_avana_Avana-2', 'avana_avana_Avana-3', 'avana_avana_Avana-4'], dtype='object')) negative control sgRNAs were found to be systematically over- or under-represented in the screens and excluded.\n",
      "  warn(\"Between %i (batch=%r) and %i (batch=%r) negative control sgRNAs were found to be \\\n",
      "2024-06-26 13:56:32.379048: W tensorflow/c/c_api.cc:304] Operation '{name:'excess_variance_1/avana/Assign' id:3132 op device:{requested: '', assigned: ''} def:{{{node excess_variance_1/avana/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](excess_variance_1/avana, excess_variance_1/avana/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "built core gene effect: 44 cell lines by 883 genes\n",
      "\n",
      "\n",
      "Connecting graph nodes in model\n",
      "\n",
      "building effective days\n",
      "built effective days, shapes {'avana': [92, 883]}\n",
      "\n",
      "building gene effect growth graph nodes\n",
      "built gene effect growth graph nodes, shapes {'avana': [44, 883]}\n",
      "\n",
      "building combined efficacy\n",
      "built combined efficacy, shape {'avana': [44, 3474]}\n",
      "built expanded combined efficacy, shapes {'avana': [92, 3474]}\n",
      "\n",
      "building growth estimates of edited cells and overall estimates of fold change in guide abundance\n",
      "built growth and change\n",
      "\n",
      "building unnormalized estimates of final abundance\n",
      "built unnormalized abundance\n",
      "\n",
      "building normalized estimates of final abundance\n",
      "built normalized abundance\n",
      "initializing precost variables\n",
      "\n",
      "\n",
      "Building all costs\n",
      "\n",
      "assembling guide efficacy regularization\n",
      "building smoothed regularization\n",
      "\n",
      "building t0 reads regularization/cost\n",
      "\n",
      "building NB2 cost\n",
      "building other regularizations\n",
      "\n",
      "Creating optimizer\n",
      "initializing rest of graph\n",
      "\tverifying graph integrity\n",
      "verifying user inputs\n",
      "verifying variables\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:56:32.822593: W tensorflow/c/c_api.cc:304] Operation '{name:'GE_1/library_effect/avana/Assign' id:3275 op device:{requested: '', assigned: ''} def:{{{node GE_1/library_effect/avana/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](GE_1/library_effect/avana, GE_1/library_effect/avana/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "verifying calculated terms\n",
      "\tavana _gene_effect\n",
      "\tavana _selected_efficacies\n",
      "\tavana_predicted_readcounts_unscaled\n",
      "\tavana _predicted_readcounts\n",
      "\tavana _normalized_readcounts\n",
      "\tavana _cost_presum\n",
      "sess run\n",
      "\tavana _cost\n",
      "\tavana _full_costs\n",
      "ready to train\n",
      "assigning trained parameters\n",
      "\tlibrary effect\n",
      "\tgene effect\n",
      "\tguide efficacy\n",
      "\tcell efficacy\n",
      "\tcell growth rate\n",
      "\tscreen excess variance\n",
      "\tscreen delay\n",
      "\tt0 offset\n",
      "Complete.\n",
      "Cost when saved: 0.212193, cost now: 0.212189\n",
      "Full cost when saved: 0.236633, full cost now: 0.236632\n"
     ]
    }
   ],
   "source": [
    "model_restored = chronos.load_saved_model(savedir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "trained model cost: 0.212193\n",
      "restored model cost: 0.212189\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:56:33.757764: W tensorflow/c/c_api.cc:304] Operation '{name:'GE/library_effect/avana/Adam_1/Assign' id:1653 op device:{requested: '', assigned: ''} def:{{{node GE/library_effect/avana/Adam_1/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](GE/library_effect/avana/Adam_1, GE/library_effect/avana/Adam_1/Initializer/zeros)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    }
   ],
   "source": [
    "print(\"trained model cost: %f\\nrestored model cost: %f\" % (model.cost, model_restored.cost))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most important file for most use cases is gene_effect.hdf5, which holds Chronos' estimate of the relative change in growth rate caused by gene knockouts. Negative values indicate inhibitory effects. You can also access the gene effect (and other parameters) from the trained model directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>gene</th>\n",
       "      <th>A1CF (29974)</th>\n",
       "      <th>A2M (2)</th>\n",
       "      <th>A2ML1 (144568)</th>\n",
       "      <th>A3GALT2 (127550)</th>\n",
       "      <th>A4GALT (53947)</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cell_line_name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ACH-000004</th>\n",
       "      <td>0.339446</td>\n",
       "      <td>0.081040</td>\n",
       "      <td>0.282277</td>\n",
       "      <td>0.464840</td>\n",
       "      <td>0.220532</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACH-000074</th>\n",
       "      <td>0.400802</td>\n",
       "      <td>0.479782</td>\n",
       "      <td>0.665330</td>\n",
       "      <td>0.390372</td>\n",
       "      <td>0.323564</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACH-000147</th>\n",
       "      <td>0.392906</td>\n",
       "      <td>0.332140</td>\n",
       "      <td>0.446861</td>\n",
       "      <td>0.144996</td>\n",
       "      <td>0.436670</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACH-000168</th>\n",
       "      <td>0.359098</td>\n",
       "      <td>0.387224</td>\n",
       "      <td>0.545619</td>\n",
       "      <td>0.307171</td>\n",
       "      <td>0.278241</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "gene            A1CF (29974)   A2M (2)  A2ML1 (144568)  A3GALT2 (127550)  \\\n",
       "cell_line_name                                                             \n",
       "ACH-000004          0.339446  0.081040        0.282277          0.464840   \n",
       "ACH-000074          0.400802  0.479782        0.665330          0.390372   \n",
       "ACH-000147          0.392906  0.332140        0.446861          0.144996   \n",
       "ACH-000168          0.359098  0.387224        0.545619          0.307171   \n",
       "\n",
       "gene            A4GALT (53947)  \n",
       "cell_line_name                  \n",
       "ACH-000004            0.220532  \n",
       "ACH-000074            0.323564  \n",
       "ACH-000147            0.436670  \n",
       "ACH-000168            0.278241  "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gene_effects = model.gene_effect\n",
    "\n",
    "gene_effects.iloc[:4, :5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copy Number Correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have gene-level copy number calls, Chronos includes an option to correct gene effect scores after the fact. This works best if the data has been scaled, as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A1CF (29974)</th>\n",
       "      <th>A2M (2)</th>\n",
       "      <th>A2ML1 (144568)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ACH-001636</th>\n",
       "      <td>1.110663</td>\n",
       "      <td>1.047670</td>\n",
       "      <td>1.047670</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACH-000784</th>\n",
       "      <td>1.193826</td>\n",
       "      <td>1.161938</td>\n",
       "      <td>1.109943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACH-000147</th>\n",
       "      <td>1.365585</td>\n",
       "      <td>0.499068</td>\n",
       "      <td>0.499068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ACH-000657</th>\n",
       "      <td>1.033547</td>\n",
       "      <td>1.040547</td>\n",
       "      <td>1.040547</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            A1CF (29974)   A2M (2)  A2ML1 (144568)\n",
       "ACH-001636      1.110663  1.047670        1.047670\n",
       "ACH-000784      1.193826  1.161938        1.109943\n",
       "ACH-000147      1.365585  0.499068        0.499068\n",
       "ACH-000657      1.033547  1.040547        1.040547"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cn = chronos.read_hdf5(\"Data/SampleData/OmicsCNGene.hdf5\")\n",
    "cn.iloc[:4, :3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, we don't have copy number calls for one of the genes targeted by the Avana library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing 1 genes from gene_effect in copy_number.\n",
      "Examples: ['POU2AF3 (120376)']\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    corrected, shifts = chronos.alternate_CN(gene_effects, cn)\n",
    "except ValueError as e:\n",
    "    print(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could choose to drop these genes. Instead, we'll assume normal ploidy (=1, in the current CCLE convention) for them and fill in the CN matrix accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in set(gene_effects.columns) - set(cn.columns):\n",
    "    cn[col] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Fitting cell line group 1 of 1\n",
      "finding low CN gene effect shifts\n",
      "smoothing and interpolating cutting toxicity for all genes\n",
      "constructed spline matrix of shape 38852, 105\n",
      "\tcost: 0.041780403256920576\n",
      "\tcost: 0.04083605505101598\n",
      "\tcost: 0.04069242969101133\n",
      "\tcost: 0.040604737769590823\n",
      "\tcost: 0.04054353723444347\n",
      "\tcost: 0.04049997688898529\n",
      "generating matrix\n"
     ]
    }
   ],
   "source": [
    "corrected, shifts = chronos.alternate_CN(gene_effects, cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `shifts` dataframe contains some information about the inferred CN effect, while `corrected` contains the corrected gene effects matrix. Overall, gene effect matrices will change little after correction, since most genes in most lines are near diploid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll write the corrected dataframe to the saved directory we made earlier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "chronos.write_hdf5(corrected, os.path.join(savedir, \"gene_effect_corrected.hdf5\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QC report"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `dataset_qc_report` in the `reports` module of Chronos presents a variety of QC metrics and interrogates some specific examples. The report minimally requires a set of positive and negative control genes. To get the full report requires copy number, mutation data, expression data, a list of expression addictions (genes which are dependencies in highly expressing lines), and oncogenic mutations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we'll load an annotated DepMap MAF file (subsetted to our cell lines). We'll select gain of function cancer driver events from it and generate a binary mutation matrix. We have a prior belief that cell lines with driver gain of function mutation events will be dependent on the mutated gene, so this matrix will be used by the QC report to assess our ability to identify selective dependencies. Specifically, we expect the oncogenes in this matrix to be dependencies in cell lines where the matrix is `True`, and not otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "maf = pd.read_csv(\"Data/SampleData/OmicsSomaticMutations.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "cancer_relevant = maf[\n",
    "  (\n",
    "      maf.Driver | maf.LikelyDriver  \n",
    "  ) & (\n",
    "      maf.LikelyGoF\n",
    "  )\n",
    "]\n",
    "\n",
    "\n",
    "cancer_relevant = cancer_relevant[~cancer_relevant.duplicated(subset=[\"ModelID\", \"Gene\"])]\n",
    "\n",
    "cancer_relevant['truecol'] = True\n",
    "\n",
    "gof_matrix_base = pd.pivot(cancer_relevant, index=\"ModelID\", columns=\"Gene\", values=\"truecol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to evaluate selective dependencies is using expression addictions, a common pattern in which a gene is a stronger dependency in lines with higher expression. We'll use a list derived from DepMap RNAi (Tsherniak et al., Cell 2017), and subset our expression matrix to match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "expression_addictions = pd.read_csv(\"Data/SampleData/RNAiExpressionAddictions.csv\")['Gene']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "addiction_expressions = chronos.read_hdf5(\"Data/SampleData/OmicsExpressionProteinCodingGenesTPMLogp1.hdf5\")[\n",
    "    expression_addictions\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we're ready to run the QC report on Chronos' results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading data from Data/Achilles_run\n",
      "plotting global control separation\n",
      "plotting selective dependency separation\n",
      "plotting gene effect mean relationships\n",
      "plotting copy number effect\n",
      "plotting screen efficacy and growth rate\n",
      "plotting readcount predictions\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/reports.py:436: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plotting LFC predictions\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/reports.py:436: UserWarning: The figure layout has changed to tight\n",
      "  plt.tight_layout()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "plotting difference from naive gene score\n",
      "summarizing\n",
      "plotting genes with low agreement with naive gene effect\n",
      "\tTNPO3 (23534)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for TNPO3 (23534), ACH-000263:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCAGTTGTTACAGATCCGGC    avGuide1\n",
      "GACTGCCTATCATATGGCCG    avGuide2\n",
      "TACCACGACCCAGATCCCAG    avGuide3\n",
      "TTACATGGTCTGGTTCCAGC    avGuide4\n",
      "dtype: object\n",
      "KASUMI-1-311cas9_RepA_p6_Avana-3    avRep1\n",
      "KASUMI-1-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for TNPO3 (23534), ACH-000527:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCAGTTGTTACAGATCCGGC    avGuide1\n",
      "GACTGCCTATCATATGGCCG    avGuide2\n",
      "TACCACGACCCAGATCCCAG    avGuide3\n",
      "TTACATGGTCTGGTTCCAGC    avGuide4\n",
      "dtype: object\n",
      "OVISE-311cas9 Rep A p6_Avana-2    avRep1\n",
      "OVISE-311cas9 Rep B p6_Avana-2    avRep2\n",
      "OVISE-311cas9 Rep C p6_Avana-2    avRep3\n",
      "dtype: object\n",
      "\tNIP7 (51388)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for NIP7 (51388), ACH-001303:\n",
      "avana    av\n",
      "dtype: object\n",
      "CTTGTCCCCGGAAATATTGG    avGuide1\n",
      "GCTGGCCGCCAATATTTCCG    avGuide2\n",
      "GTTTCCGTCTGCACAACGAC    avGuide3\n",
      "TATCTTCTCAAACATGACAC    avGuide4\n",
      "dtype: object\n",
      "NB1643-311Cas9_RepA_p4_Avana-3    avRep1\n",
      "NB1643-311Cas9_RepB_p4_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for NIP7 (51388), ACH-000912:\n",
      "avana    av\n",
      "dtype: object\n",
      "CTTGTCCCCGGAAATATTGG    avGuide1\n",
      "GCTGGCCGCCAATATTTCCG    avGuide2\n",
      "GTTTCCGTCTGCACAACGAC    avGuide3\n",
      "TATCTTCTCAAACATGACAC    avGuide4\n",
      "dtype: object\n",
      "NCI-H2286-311caa9_RepA_p6_Avana-3    avRep1\n",
      "NCI-H2286-311caa9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tESR1 (2099)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for ESR1 (2099), ACH-000147:\n",
      "avana    av\n",
      "dtype: object\n",
      "CGCCGTGTACAACTACCCCG    avGuide1\n",
      "GCACCATTGATAAAAACAGG    avGuide2\n",
      "GTAGACCTGCGCGTTGGCGG    avGuide3\n",
      "TCGCCTTTCCTGCAGCCCCA    avGuide4\n",
      "dtype: object\n",
      "T47D-311Cas9-RepA-p6_Avana-4    avRep1\n",
      "T47D-311Cas9-RepB-p6_Avana-4    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for ESR1 (2099), ACH-001842:\n",
      "avana    av\n",
      "dtype: object\n",
      "CGCCGTGTACAACTACCCCG    avGuide1\n",
      "GCACCATTGATAAAAACAGG    avGuide2\n",
      "GTAGACCTGCGCGTTGGCGG    avGuide3\n",
      "TCGCCTTTCCTGCAGCCCCA    avGuide4\n",
      "dtype: object\n",
      "ICC2-311cas9_RepA_p6_Avana-3    avRep1\n",
      "ICC2-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tPRMT5 (10419)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for PRMT5 (10419), ACH-000988:\n",
      "avana    av\n",
      "dtype: object\n",
      "ATGAACTCCCTCTTGAAACG    avGuide1\n",
      "CTAGTCATCCCGGAGAAGCA    avGuide2\n",
      "GGAGAAAAACCCAAATGCCG    avGuide3\n",
      "GGTTGCTACTCACGTCACCA    avGuide4\n",
      "dtype: object\n",
      "MFE-319-311Cas9_RepA_p6_Avana-3    avRep1\n",
      "MFE-319-311Cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for PRMT5 (10419), ACH-000935:\n",
      "avana    av\n",
      "dtype: object\n",
      "ATGAACTCCCTCTTGAAACG    avGuide1\n",
      "CTAGTCATCCCGGAGAAGCA    avGuide2\n",
      "GGAGAAAAACCCAAATGCCG    avGuide3\n",
      "GGTTGCTACTCACGTCACCA    avGuide4\n",
      "dtype: object\n",
      "MDST8-311Cas9 Rep A p6_Avana-3    avRep1\n",
      "MDST8-311Cas9 Rep C p6_Avana-3    avRep2\n",
      "MDST8-311Cas9 Rep D p6_Avana-3    avRep3\n",
      "dtype: object\n",
      "\tJAK3 (3718)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for JAK3 (3718), ACH-000406:\n",
      "avana    av\n",
      "dtype: object\n",
      "GGAAGCTGTCAAAGTCCTGG    avGuide1\n",
      "GGAGCGGCTGGATCCAGCCG    avGuide2\n",
      "TCCTTCACCAAGATTTACCG    avGuide3\n",
      "TGCCGCCGTCACCAGCCACG    avGuide4\n",
      "dtype: object\n",
      "U937-101Cas9 Rep A p6_Avana-3    avRep1\n",
      "U937-101Cas9 Rep C p6_Avana-3    avRep2\n",
      "U937-101Cas9 Rep D p6_Avana-3    avRep3\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for JAK3 (3718), ACH-000479:\n",
      "avana    av\n",
      "dtype: object\n",
      "GGAAGCTGTCAAAGTCCTGG    avGuide1\n",
      "GGAGCGGCTGGATCCAGCCG    avGuide2\n",
      "TCCTTCACCAAGATTTACCG    avGuide3\n",
      "TGCCGCCGTCACCAGCCACG    avGuide4\n",
      "dtype: object\n",
      "KNS-81-311cas9_RepA_p5_Avana-3    avRep1\n",
      "KNS-81-311cas9_RepB_p5_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tJAK2 (3717)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for JAK2 (3717), ACH-000004:\n",
      "avana    av\n",
      "dtype: object\n",
      "ATAATCATACCTGAGCAAAG    avGuide1\n",
      "ATCTGCCTCAGATTTCCCAA    avGuide2\n",
      "CTGCCACTGCAATACCAACG    avGuide3\n",
      "GAGGTGCTACTTCTTTACAG    avGuide4\n",
      "dtype: object\n",
      "HEL-311Cas9_RepA_p4_Avana-3    avRep1\n",
      "HEL-311Cas9_RepB_p4_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for JAK2 (3717), ACH-000396:\n",
      "avana    av\n",
      "dtype: object\n",
      "ATAATCATACCTGAGCAAAG    avGuide1\n",
      "ATCTGCCTCAGATTTCCCAA    avGuide2\n",
      "CTGCCACTGCAATACCAACG    avGuide3\n",
      "GAGGTGCTACTTCTTTACAG    avGuide4\n",
      "dtype: object\n",
      "J82-311Cas9_RepA_p5_Avana-3    avRep1\n",
      "J82-311Cas9_RepB_p5_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tSSU72 (29101)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for SSU72 (29101), ACH-001736:\n",
      "avana    av\n",
      "dtype: object\n",
      "GAAATAAGAGAATCAAGCCC    avGuide1\n",
      "GAGTCCCTGTTCCAAAGGAT    avGuide2\n",
      "GCGTCCGATCCTTTGGAACA    avGuide3\n",
      "TTCTGGAATCTTTCTGGCCG    avGuide4\n",
      "dtype: object\n",
      "HB11;19-311CAS9_RepA_p6_Avana-3    avRep1\n",
      "HB11;19-311CAS9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for SSU72 (29101), ACH-000168:\n",
      "avana    av\n",
      "dtype: object\n",
      "GAAATAAGAGAATCAAGCCC    avGuide1\n",
      "GAGTCCCTGTTCCAAAGGAT    avGuide2\n",
      "GCGTCCGATCCTTTGGAACA    avGuide3\n",
      "TTCTGGAATCTTTCTGGCCG    avGuide4\n",
      "dtype: object\n",
      "NOMO-1-311Cas9_RepA_p4_Avana-2    avRep1\n",
      "NOMO-1-311Cas9_RepB_p4_Avana-2    avRep2\n",
      "dtype: object\n",
      "\tMCM6 (4175)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for MCM6 (4175), ACH-000788:\n",
      "avana    av\n",
      "dtype: object\n",
      "GGCAATGATGAAGTAAAACG    avGuide1\n",
      "GTTCACCCAGAGCTTGTGAG    avGuide2\n",
      "TCAAGGCCCGACACAGGTAA    avGuide3\n",
      "TGTTGCGCCAACCAACCCAA    avGuide4\n",
      "dtype: object\n",
      "A2058-311cas9_RepA_p6_Avana-3    avRep1\n",
      "A2058-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for MCM6 (4175), ACH-000750:\n",
      "avana    av\n",
      "dtype: object\n",
      "GGCAATGATGAAGTAAAACG    avGuide1\n",
      "GTTCACCCAGAGCTTGTGAG    avGuide2\n",
      "TCAAGGCCCGACACAGGTAA    avGuide3\n",
      "TGTTGCGCCAACCAACCCAA    avGuide4\n",
      "dtype: object\n",
      "LOXIMVI-311Cas9_RepA_p6_Avana-2    avRep1\n",
      "LOXIMVI-311Cas9_RepB_p6_Avana-2    avRep2\n",
      "dtype: object\n",
      "\tNEDD1 (121441)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for NEDD1 (121441), ACH-000614:\n",
      "avana    av\n",
      "dtype: object\n",
      "GGCTGACACTCCTCTAACTG    avGuide1\n",
      "TAGAACTGTGGCAATGGACG    avGuide2\n",
      "TCTATGTATTTGGTAAGCGG    avGuide3\n",
      "TGGCTATTGGATCTTCCCGG    avGuide4\n",
      "dtype: object\n",
      "RVH421-311Cas9_RepA_p5_Avana-3        avRep1\n",
      "RVH421-311Cas9_RepB_p5_Avana-3        avRep2\n",
      "RVH421RPMI-311Cas9-RepA-p6_Avana-4    avRep3\n",
      "RVH421RPMI-311Cas9-RepB-p6_Avana-4    avRep4\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for NEDD1 (121441), ACH-001843:\n",
      "avana    av\n",
      "dtype: object\n",
      "GGCTGACACTCCTCTAACTG    avGuide1\n",
      "TAGAACTGTGGCAATGGACG    avGuide2\n",
      "TCTATGTATTTGGTAAGCGG    avGuide3\n",
      "TGGCTATTGGATCTTCCCGG    avGuide4\n",
      "dtype: object\n",
      "ICC3-311cas9_RepA_p6_Avana-3    avRep1\n",
      "ICC3-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tRFT1 (91869)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for RFT1 (91869), ACH-001674:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCCCACTGAGACATGCTCTG    avGuide1\n",
      "GCTTGTGCCAAGACCCAAAA    avGuide2\n",
      "GGAGGCTATGGGCAGCCAGG    avGuide3\n",
      "GGTGATTGCAGAGAGCCTGT    avGuide4\n",
      "dtype: object\n",
      "TGW-311Cas9-RepA-P6_Avana-4    avRep1\n",
      "TGW-311Cas9-RepB-P6_Avana-4    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for RFT1 (91869), ACH-001636:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCCCACTGAGACATGCTCTG    avGuide1\n",
      "GCTTGTGCCAAGACCCAAAA    avGuide2\n",
      "GGAGGCTATGGGCAGCCAGG    avGuide3\n",
      "GGTGATTGCAGAGAGCCTGT    avGuide4\n",
      "dtype: object\n",
      "RAMOS-311cas9-RepA-p6_Avana-4    avRep1\n",
      "RAMOS-311cas9-RepB-p6_Avana-4    avRep2\n",
      "dtype: object\n",
      "\tRRP1 (8568)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for RRP1 (8568), ACH-001736:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCTTCTGCAGAATTGCTGCC    avGuide1\n",
      "CTTCATCGAGATCTTCCTGG    avGuide2\n",
      "GCAGGTGACCCGGGACCGGG    avGuide3\n",
      "GCGCCTGGATAAATTCTACA    avGuide4\n",
      "dtype: object\n",
      "HB11;19-311CAS9_RepA_p6_Avana-3    avRep1\n",
      "HB11;19-311CAS9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for RRP1 (8568), ACH-000750:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCTTCTGCAGAATTGCTGCC    avGuide1\n",
      "CTTCATCGAGATCTTCCTGG    avGuide2\n",
      "GCAGGTGACCCGGGACCGGG    avGuide3\n",
      "GCGCCTGGATAAATTCTACA    avGuide4\n",
      "dtype: object\n",
      "LOXIMVI-311Cas9_RepA_p6_Avana-2    avRep1\n",
      "LOXIMVI-311Cas9_RepB_p6_Avana-2    avRep2\n",
      "dtype: object\n",
      "\tSEC61A1 (29927)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for SEC61A1 (29927), ACH-000458:\n",
      "avana    av\n",
      "dtype: object\n",
      "ACTACTGTCAACACTGGCCG    avGuide1\n",
      "CTGGCTTCTGAATTTCCGGC    avGuide2\n",
      "GCACACTAAGAAGATAAAGA    avGuide3\n",
      "GTGATTCTAGCCTCTAACAG    avGuide4\n",
      "dtype: object\n",
      "CJM-311Cas9_RepA_p7_Avana-3    avRep1\n",
      "CJM-311Cas9_RepB_p7_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for SEC61A1 (29927), ACH-000147:\n",
      "avana    av\n",
      "dtype: object\n",
      "ACTACTGTCAACACTGGCCG    avGuide1\n",
      "CTGGCTTCTGAATTTCCGGC    avGuide2\n",
      "GCACACTAAGAAGATAAAGA    avGuide3\n",
      "GTGATTCTAGCCTCTAACAG    avGuide4\n",
      "dtype: object\n",
      "T47D-311Cas9-RepA-p6_Avana-4    avRep1\n",
      "T47D-311Cas9-RepB-p6_Avana-4    avRep2\n",
      "dtype: object\n",
      "\tCLEC2A (387836)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for CLEC2A (387836), ACH-000004:\n",
      "avana    av\n",
      "dtype: object\n",
      "ATTAATCCAGAGCTGCGGGA    avGuide1\n",
      "CCATGATTAATCCAGAGCTG    avGuide2\n",
      "GAGCTGATGGCTTCATACAT    avGuide3\n",
      "GATGGCTTCATACATCGGAT    avGuide4\n",
      "dtype: object\n",
      "HEL-311Cas9_RepA_p4_Avana-3    avRep1\n",
      "HEL-311Cas9_RepB_p4_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for CLEC2A (387836), ACH-001636:\n",
      "avana    av\n",
      "dtype: object\n",
      "ATTAATCCAGAGCTGCGGGA    avGuide1\n",
      "CCATGATTAATCCAGAGCTG    avGuide2\n",
      "GAGCTGATGGCTTCATACAT    avGuide3\n",
      "GATGGCTTCATACATCGGAT    avGuide4\n",
      "dtype: object\n",
      "RAMOS-311cas9-RepA-p6_Avana-4    avRep1\n",
      "RAMOS-311cas9-RepB-p6_Avana-4    avRep2\n",
      "dtype: object\n",
      "\tGUCA2A (2980)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for GUCA2A (2980), ACH-000996:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCAGGAGCCCCAGGAGCCCA    avGuide1\n",
      "GCACAGTGCGGAGAGCAGGA    avGuide2\n",
      "GGCCCCAAGGAGGCACAGTG    avGuide3\n",
      "TGGGCTCCTGGGGCTCCTGG    avGuide4\n",
      "dtype: object\n",
      "HEC-251-311Cas9_RepA_p5_Avana-3    avRep1\n",
      "HEC-251-311Cas9_RepB_p5_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for GUCA2A (2980), ACH-001636:\n",
      "avana    av\n",
      "dtype: object\n",
      "CCAGGAGCCCCAGGAGCCCA    avGuide1\n",
      "GCACAGTGCGGAGAGCAGGA    avGuide2\n",
      "GGCCCCAAGGAGGCACAGTG    avGuide3\n",
      "TGGGCTCCTGGGGCTCCTGG    avGuide4\n",
      "dtype: object\n",
      "RAMOS-311cas9-RepA-p6_Avana-4    avRep1\n",
      "RAMOS-311cas9-RepB-p6_Avana-4    avRep2\n",
      "dtype: object\n",
      "\tBTK (695)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for BTK (695), ACH-001533:\n",
      "avana    av\n",
      "dtype: object\n",
      "AATGATCTACAGCTGCGGAA    avGuide1\n",
      "ATAGTAGGAGAGTTTGTGCA    avGuide2\n",
      "CTGTGTTTGCTAAATCCACA    avGuide3\n",
      "GATGCTCTCCAGAATCACTG    avGuide4\n",
      "dtype: object\n",
      "KARPAS1718-311cas9-RepA-p6_Avana-4    avRep1\n",
      "KARPAS1718-311cas9-RepB-p6_Avana-4    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for BTK (695), ACH-000183:\n",
      "avana    av\n",
      "dtype: object\n",
      "AATGATCTACAGCTGCGGAA    avGuide1\n",
      "ATAGTAGGAGAGTTTGTGCA    avGuide2\n",
      "CTGTGTTTGCTAAATCCACA    avGuide3\n",
      "GATGCTCTCCAGAATCACTG    avGuide4\n",
      "dtype: object\n",
      "L363-311Cas9_RepA_p6_Avana-3    avRep1\n",
      "L363-311Cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tPOU2AF3 (120376)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for POU2AF3 (120376), ACH-001843:\n",
      "avana    av\n",
      "dtype: object\n",
      "CTTTGTGTGGAATTGAGCAA    avGuide1\n",
      "GCTCTGGTTAAAGCAGAAAG    avGuide2\n",
      "TCTGGTGTCCAGAGAAGCAG    avGuide3\n",
      "TTCTGGTGTCCAGAGAAGCA    avGuide4\n",
      "dtype: object\n",
      "ICC3-311cas9_RepA_p6_Avana-3    avRep1\n",
      "ICC3-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for POU2AF3 (120376), ACH-000750:\n",
      "avana    av\n",
      "dtype: object\n",
      "CTTTGTGTGGAATTGAGCAA    avGuide1\n",
      "GCTCTGGTTAAAGCAGAAAG    avGuide2\n",
      "TCTGGTGTCCAGAGAAGCAG    avGuide3\n",
      "TTCTGGTGTCCAGAGAAGCA    avGuide4\n",
      "dtype: object\n",
      "LOXIMVI-311Cas9_RepA_p6_Avana-2    avRep1\n",
      "LOXIMVI-311Cas9_RepB_p6_Avana-2    avRep2\n",
      "dtype: object\n",
      "\tPRDM1 (639)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for PRDM1 (639), ACH-000183:\n",
      "avana    av\n",
      "dtype: object\n",
      "CTTGAACGCGTCCTACGGCA    avGuide1\n",
      "GGATGGGGTAAACGACCCGA    avGuide2\n",
      "TTACAATTCATGCCGTAGGG    avGuide3\n",
      "TTGGACAGATCTATTCCAGA    avGuide4\n",
      "dtype: object\n",
      "L363-311Cas9_RepA_p6_Avana-3    avRep1\n",
      "L363-311Cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for PRDM1 (639), ACH-000479:\n",
      "avana    av\n",
      "dtype: object\n",
      "CTTGAACGCGTCCTACGGCA    avGuide1\n",
      "GGATGGGGTAAACGACCCGA    avGuide2\n",
      "TTACAATTCATGCCGTAGGG    avGuide3\n",
      "TTGGACAGATCTATTCCAGA    avGuide4\n",
      "dtype: object\n",
      "KNS-81-311cas9_RepA_p5_Avana-3    avRep1\n",
      "KNS-81-311cas9_RepB_p5_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tFOXR1 (283150)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for FOXR1 (283150), ACH-000946:\n",
      "avana    av\n",
      "dtype: object\n",
      "GAGACCTCCAGCTTTCCAGG    avGuide1\n",
      "GGAAGATGCCAGCTGCTCAG    avGuide2\n",
      "TGAGACCTCCAGCTTTCCAG    avGuide3\n",
      "TGGGATTTACCCACATCCAG    avGuide4\n",
      "dtype: object\n",
      "HEC265-311Cas9_RepA_p6_Avana-3    avRep1\n",
      "HEC265-311Cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for FOXR1 (283150), ACH-000750:\n",
      "avana    av\n",
      "dtype: object\n",
      "GAGACCTCCAGCTTTCCAGG    avGuide1\n",
      "GGAAGATGCCAGCTGCTCAG    avGuide2\n",
      "TGAGACCTCCAGCTTTCCAG    avGuide3\n",
      "TGGGATTTACCCACATCCAG    avGuide4\n",
      "dtype: object\n",
      "LOXIMVI-311Cas9_RepA_p6_Avana-2    avRep1\n",
      "LOXIMVI-311Cas9_RepB_p6_Avana-2    avRep2\n",
      "dtype: object\n",
      "\tPFDN2 (5202)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for PFDN2 (5202), ACH-000988:\n",
      "avana    av\n",
      "dtype: object\n",
      "AACTCGTAAGTGCTACCGCA    avGuide1\n",
      "CCAGCACTCCTCCAACCATG    avGuide2\n",
      "GCGAAGATGGCGGAGAACAG    avGuide3\n",
      "GGTCGCGCCGGCAAGAGCAG    avGuide4\n",
      "dtype: object\n",
      "MFE-319-311Cas9_RepA_p6_Avana-3    avRep1\n",
      "MFE-319-311Cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for PFDN2 (5202), ACH-000334:\n",
      "avana    av\n",
      "dtype: object\n",
      "AACTCGTAAGTGCTACCGCA    avGuide1\n",
      "CCAGCACTCCTCCAACCATG    avGuide2\n",
      "GCGAAGATGGCGGAGAACAG    avGuide3\n",
      "GGTCGCGCCGGCAAGAGCAG    avGuide4\n",
      "dtype: object\n",
      "DB-311Cas9_RepA_p3_Avana-3    avRep1\n",
      "DB-311Cas9_RepB_p3_Avana-3    avRep2\n",
      "dtype: object\n",
      "\tLRIT2 (340745)\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for LRIT2 (340745), ACH-001843:\n",
      "avana    av\n",
      "dtype: object\n",
      "GCCAGGATGCACTTGCTCAG    avGuide1\n",
      "GGCCTGTGCCAAGCATCGCA    avGuide2\n",
      "GGCGACAGTCACATACCCAG    avGuide3\n",
      "GGTATGTGACTGTCGCCTAA    avGuide4\n",
      "dtype: object\n",
      "ICC3-311cas9_RepA_p6_Avana-3    avRep1\n",
      "ICC3-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "avana    avana\n",
      "dtype: object avana    avana\n",
      "dtype: object\n",
      "Guide and replicate key for LRIT2 (340745), ACH-001129:\n",
      "avana    av\n",
      "dtype: object\n",
      "GCCAGGATGCACTTGCTCAG    avGuide1\n",
      "GGCCTGTGCCAAGCATCGCA    avGuide2\n",
      "GGCGACAGTCACATACCCAG    avGuide3\n",
      "GGTATGTGACTGTCGCCTAA    avGuide4\n",
      "dtype: object\n",
      "MONO-MAC1-311cas9_RepA_p6_Avana-3    avRep1\n",
      "MONO-MAC1-311cas9_RepB_p6_Avana-3    avRep2\n",
      "dtype: object\n",
      "building report\n"
     ]
    }
   ],
   "source": [
    "metrics = reports.dataset_qc_report(\"ChronosAvana\", savedir, \n",
    "                          common_essentials, nonessentials,\n",
    "                          gof_matrix_base, addiction_expressions,\n",
    "                          cn, directory=\"Data/reports\",\n",
    "                          gene_effect_file=\"gene_effect_corrected.hdf5\"\n",
    "                         )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running with multiple libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add Sanger's [Project Score](https://www.nature.com/articles/s41586-019-1103-9) data (screened with the KY library) and run Chronos jointly on it and the Avana data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "ky_guide_map = pd.read_csv(\"./Data/SampleData/KYGuideMap.csv\")\n",
    "ky_sequence_map = pd.read_csv(\"./Data/SampleData/KYSequenceMap.csv\")\n",
    "ky_readcounts = chronos.read_hdf5(\"./Data/SampleData/KYReadcounts.hdf5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "ky_positive_controls = ky_guide_map.sgrna[ky_guide_map.gene.isin(common_essentials)]\n",
    "ky_negative_controls = ky_guide_map.sgrna[ky_guide_map.gene.isin(nonessentials)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the call signature of Chronos with multiple libraries is constructed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normalizing readcounts\n",
      "\n",
      "\n",
      "Finding all unique guides and genes\n",
      "found 3474 unique guides and 883 unique genes in avana\n",
      "found 4084 unique guides and 833 unique genes in ky\n",
      "found 7558 unique guides and 887 unique genes overall\n",
      "\n",
      "finding guide-gene mapping indices\n",
      "\n",
      "finding all unique sequenced replicates, cell lines, and pDNA batches\n",
      "found 92 unique sequences (excluding pDNA) and 44 unique cell lines in avana\n",
      "found 63 unique sequences (excluding pDNA) and 23 unique cell lines in ky\n",
      "found 155 unique replicates and 58 unique cell lines overall\n",
      "\n",
      "finding replicate-cell line mappings indices\n",
      "\n",
      "finding replicate-pDNA mappings indices\n",
      "\n",
      "\n",
      "assigning float constants\n",
      "Estimating or aligning variances\n",
      "\tEstimating excess variance (alpha) for avana\n",
      "\tEstimating excess variance (alpha) for ky\n",
      "Creating excess variance tensors\n",
      "\tCreated excess variance tensor for avana with shape [92, 1]\n",
      "\tCreated excess variance tensor for ky with shape [63, 1]\n",
      "initializing graph\n",
      "\n",
      "building gene effect mask\n",
      "\n",
      "building doubling vectors\n",
      "made days vector of shape [92, 1] for avana\n",
      "made days vector of shape [63, 1] for ky\n",
      "\n",
      "building late observed timepoints\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/model.py:377: UserWarning: Between 0 (batch=Index(['avana_Avana-2', 'avana_Avana-3', 'avana_Avana-4'], dtype='object')) and 0 (batch=Index(['avana_Avana-2', 'avana_Avana-3', 'avana_Avana-4'], dtype='object')) negative control sgRNAs were found to be systematically over- or under-represented in the screens and excluded.\n",
      "  warn(\"Between %i (batch=%r) and %i (batch=%r) negative control sgRNAs were found to be \\\n",
      "/Users/dempster/packages/chronos_public/chronos/model.py:377: UserWarning: Between 4 (batch=Index(['ky_KY-1'], dtype='object')) and 39 (batch=Index(['ky_KY-2'], dtype='object')) negative control sgRNAs were found to be systematically over- or under-represented in the screens and excluded.\n",
      "  warn(\"Between %i (batch=%r) and %i (batch=%r) negative control sgRNAs were found to be \\\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\tbuilt normalized timepoints for avana with shape (92, 3474) (replicates X guides)\n",
      "\tbuilt normalized timepoints for ky with shape (63, 4084) (replicates X guides)\n",
      "\n",
      "building t0 reads\n",
      "\n",
      "\n",
      "Building variables\n",
      "\n",
      "building t0 reads estimate\n",
      "made t0 batch with shape (2, 4084) for avana\n",
      "made t0 batch with shape (2, 4084) for ky\n",
      "building guide efficacy\n",
      "built guide efficacy: shape [1, 7559]\n",
      "building growth rate\n",
      "built growth rate: shape {'avana': [58, 1], 'ky': [58, 1]}\n",
      "\n",
      "building cell line efficacy\n",
      "built cell line efficacy: shapes {'avana': [58, 1], 'ky': [58, 1]}\n",
      "building screen delay\n",
      "built screen delay\n",
      "building gene effect\n",
      "built core gene effect: 58 cell lines by 887 genes\n",
      "\n",
      "\n",
      "Connecting graph nodes in model\n",
      "\n",
      "building effective days\n",
      "built effective days, shapes {'avana': [92, 887], 'ky': [63, 887]}\n",
      "\n",
      "building gene effect growth graph nodes\n",
      "built gene effect growth graph nodes, shapes {'avana': [58, 887], 'ky': [58, 887]}\n",
      "\n",
      "building combined efficacy\n",
      "built combined efficacy, shape {'avana': [58, 3474], 'ky': [58, 4084]}\n",
      "built expanded combined efficacy, shapes {'avana': [92, 3474], 'ky': [63, 4084]}\n",
      "\n",
      "building growth estimates of edited cells and overall estimates of fold change in guide abundance\n",
      "built growth and change\n",
      "\n",
      "building unnormalized estimates of final abundance\n",
      "built unnormalized abundance\n",
      "\n",
      "building normalized estimates of final abundance\n",
      "built normalized abundance\n",
      "initializing precost variables\n",
      "\n",
      "\n",
      "Building all costs\n",
      "\n",
      "assembling guide efficacy regularization\n",
      "building smoothed regularization\n",
      "\n",
      "building t0 reads regularization/cost\n",
      "\n",
      "building NB2 cost\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:57:29.263786: W tensorflow/c/c_api.cc:304] Operation '{name:'inferred_t0_2/base_avana/Assign' id:6377 op device:{requested: '', assigned: ''} def:{{{node inferred_t0_2/base_avana/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](inferred_t0_2/base_avana, inferred_t0_2/base_avana/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building other regularizations\n",
      "\n",
      "Creating optimizer\n",
      "initializing rest of graph\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:57:30.112467: W tensorflow/c/c_api.cc:304] Operation '{name:'GE_2/library_effect/ky/Assign' id:6592 op device:{requested: '', assigned: ''} def:{{{node GE_2/library_effect/ky/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](GE_2/library_effect/ky, GE_2/library_effect/ky/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimating initial screen efficacy and gene effect\n",
      "\t avana\n",
      "\t ky\n",
      "\tverifying graph integrity\n",
      "verifying user inputs\n",
      "verifying variables\n",
      "verifying calculated terms\n",
      "\tavana _gene_effect\n",
      "\tavana _selected_efficacies\n",
      "\tavana_predicted_readcounts_unscaled\n",
      "\tavana _predicted_readcounts\n",
      "\tavana _normalized_readcounts\n",
      "\tavana _cost_presum\n",
      "sess run\n",
      "\tavana _cost\n",
      "\tavana _full_costs\n",
      "\tky _gene_effect\n",
      "\tky _selected_efficacies\n",
      "\tky_predicted_readcounts_unscaled\n",
      "\tky _predicted_readcounts\n",
      "\tky _normalized_readcounts\n",
      "\tky _cost_presum\n",
      "sess run\n",
      "\tky _cost\n",
      "\tky _full_costs\n",
      "ready to train\n"
     ]
    }
   ],
   "source": [
    "model2 = chronos.Chronos(\n",
    "    sequence_map={\"avana\": sequence_map, 'ky': ky_sequence_map},\n",
    "    guide_gene_map={\"avana\": guide_map, 'ky': ky_guide_map},\n",
    "    readcounts={\"avana\": readcounts, 'ky': ky_readcounts},\n",
    "    negative_control_sgrnas={\"avana\": negative_controls, \"ky\": ky_negative_controls}\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NB2 cost 0.3618576638967593\n",
      "Full cost 0.5103175488433159\n",
      "relative_growth_rate\n",
      "\tavana max 1.037, min 0.96811\n",
      "\tky max 1.044, min 0.95815\n",
      "mean guide efficacy 0.9991983630066278\n",
      "t0_offset SD: [('avana', 7.040503180981407e-05), ('ky', 7.552244104967997e-05)]\n",
      "\n",
      "gene mean -0.26400131756601397\n",
      "SD of gene means 0.4294472633028505\n",
      "Mean of gene SDs 0.3446919369883686\n",
      "\n",
      "\n",
      "\n",
      "51 epochs trained, time taken 0:00:03, projected remaining 0:00:14\n",
      "NB2 cost 0.25407094118085255\n",
      "Full cost 0.3154294898357065\n",
      "relative_growth_rate\n",
      "\tavana max 1.256, min 0.53759\n",
      "\tky max 1.468, min 0.47675\n",
      "mean guide efficacy 0.9370800264295374\n",
      "t0_offset SD: [('avana', 0.14024269896206715), ('ky', 0.1386491460853011)]\n",
      "\n",
      "gene mean -0.02314857033179854\n",
      "SD of gene means 0.5424142871732163\n",
      "Mean of gene SDs 0.2769753406695022\n",
      "\n",
      "\n",
      "\n",
      "101 epochs trained, time taken 0:00:05, projected remaining 0:00:10\n",
      "NB2 cost 0.23471048316686866\n",
      "Full cost 0.25862748334718416\n",
      "relative_growth_rate\n",
      "\tavana max 1.544, min 0.52428\n",
      "\tky max 1.275, min 0.70816\n",
      "mean guide efficacy 0.8745269585011546\n",
      "t0_offset SD: [('avana', 0.13577795214869234), ('ky', 0.11442315897431432)]\n",
      "\n",
      "gene mean 0.00041532150432264687\n",
      "SD of gene means 0.6415058424029152\n",
      "Mean of gene SDs 0.21704474877209867\n",
      "\n",
      "\n",
      "\n",
      "151 epochs trained, time taken 0:00:08, projected remaining 0:00:07\n",
      "NB2 cost 0.2338909937276131\n",
      "Full cost 0.2565336029079192\n",
      "relative_growth_rate\n",
      "\tavana max 1.509, min 0.51746\n",
      "\tky max 1.383, min 0.71219\n",
      "mean guide efficacy 0.8711999362819592\n",
      "t0_offset SD: [('avana', 0.13483669730744952), ('ky', 0.11276787072494918)]\n",
      "\n",
      "gene mean -0.00016303144696584265\n",
      "SD of gene means 0.6426578400695825\n",
      "Mean of gene SDs 0.20338351270626945\n",
      "\n",
      "\n",
      "\n",
      "201 epochs trained, time taken 0:00:10, projected remaining 0:00:05\n",
      "NB2 cost 0.23347697894905325\n",
      "Full cost 0.2555518707170642\n",
      "relative_growth_rate\n",
      "\tavana max 1.502, min 0.52609\n",
      "\tky max 1.397, min 0.71241\n",
      "mean guide efficacy 0.8679009524670551\n",
      "t0_offset SD: [('avana', 0.13346522858285942), ('ky', 0.11111942424775316)]\n",
      "\n",
      "gene mean 0.0004950490623915282\n",
      "SD of gene means 0.6465781459843389\n",
      "Mean of gene SDs 0.19928061571991795\n",
      "\n",
      "\n",
      "\n",
      "251 epochs trained, time taken 0:00:12, projected remaining 0:00:02\n",
      "NB2 cost 0.23329294362318048\n",
      "Full cost 0.25520428284624064\n",
      "relative_growth_rate\n",
      "\tavana max 1.515, min 0.52302\n",
      "\tky max 1.415, min 0.69954\n",
      "mean guide efficacy 0.8650010503730644\n",
      "t0_offset SD: [('avana', 0.13292015685359096), ('ky', 0.11010597878647038)]\n",
      "\n",
      "gene mean 0.0009263036990135608\n",
      "SD of gene means 0.6473512538486412\n",
      "Mean of gene SDs 0.19798909582507038\n",
      "\n",
      "\n",
      "\n",
      "301 epochs trained, time taken 0:00:15, projected remaining 0:00:00\n",
      "NB2 cost 0.23327909878566291\n",
      "Full cost 0.2551970651785153\n",
      "relative_growth_rate\n",
      "\tavana max 1.509, min 0.51652\n",
      "\tky max 1.421, min 0.70058\n",
      "mean guide efficacy 0.863572209267812\n",
      "t0_offset SD: [('avana', 0.13329488988825747), ('ky', 0.1097737320624151)]\n",
      "\n",
      "gene mean 0.0014370125171794417\n",
      "SD of gene means 0.6474714783907901\n",
      "Mean of gene SDs 0.1977718988222438\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model2.train(301)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the gene effect now has NAs. These are cases where a cell line was only screened in one library and that library had no guides for that gene."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chronos infers library batch effects. Note that these are only inferred for genes present in all libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>library</th>\n",
       "      <th>avana</th>\n",
       "      <th>ky</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A1CF (29974)</th>\n",
       "      <td>-0.138539</td>\n",
       "      <td>0.111505</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A2M (2)</th>\n",
       "      <td>-0.282217</td>\n",
       "      <td>0.227145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A2ML1 (144568)</th>\n",
       "      <td>-0.197223</td>\n",
       "      <td>0.158737</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A3GALT2 (127550)</th>\n",
       "      <td>-0.032770</td>\n",
       "      <td>0.026375</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>A4GALT (53947)</th>\n",
       "      <td>-0.112406</td>\n",
       "      <td>0.090471</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZNF622 (90441)</th>\n",
       "      <td>-0.219404</td>\n",
       "      <td>0.176590</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZNF648 (127665)</th>\n",
       "      <td>-0.177326</td>\n",
       "      <td>0.142723</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZNF804B (219578)</th>\n",
       "      <td>-0.235083</td>\n",
       "      <td>0.189209</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZP2 (7783)</th>\n",
       "      <td>-0.238114</td>\n",
       "      <td>0.191648</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ZSWIM2 (151112)</th>\n",
       "      <td>-0.241264</td>\n",
       "      <td>0.194183</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>887 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "library              avana        ky\n",
       "gene                                \n",
       "A1CF (29974)     -0.138539  0.111505\n",
       "A2M (2)          -0.282217  0.227145\n",
       "A2ML1 (144568)   -0.197223  0.158737\n",
       "A3GALT2 (127550) -0.032770  0.026375\n",
       "A4GALT (53947)   -0.112406  0.090471\n",
       "...                    ...       ...\n",
       "ZNF622 (90441)   -0.219404  0.176590\n",
       "ZNF648 (127665)  -0.177326  0.142723\n",
       "ZNF804B (219578) -0.235083  0.189209\n",
       "ZP2 (7783)       -0.238114  0.191648\n",
       "ZSWIM2 (151112)  -0.241264  0.194183\n",
       "\n",
       "[887 rows x 2 columns]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model2.library_effect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running your screen with pretrained DepMap parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you conducted a screen in one of the DepMap integrated libraries (currently Avana, KY, or Humagne-CD), you can load parameters from the trained DepMap model and use them to process your specific screen. This gives you many of the benefits of coprocessing your screen with the complete DepMap dataset without the computational expense. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following command fetches the 22Q3 public dataset from Figshare and stores it in the Chronos package directory under Data/DepMapParameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "`url_loc` will be found relative to the chronos package directory\n",
      "'/Users/dempster/packages/chronos_public/chronos'\n",
      "Pass `relative_to_chronos=False` to make the path relative to your current working directory\n",
      "'/Users/dempster/packages/chronos_public'\n",
      "instead.\n",
      "\n",
      "`output_dir` will be found relative to the chronos package directory\n",
      "'/Users/dempster/packages/chronos_public/chronos'\n",
      "Pass `relative_to_chronos=False` to make the path relative to your current working directory\n",
      "'/Users/dempster/packages/chronos_public'\n",
      "instead.\n",
      "\n",
      "downloading files to /Users/dempster/packages/chronos_public/chronos/../Data/DepMapParameters/\n",
      "Skipping gene_effect.csv as it already exists, pass `overwrite=True` to overwrite\n",
      "Skipping guide_efficacy.csv as it already exists, pass `overwrite=True` to overwrite\n",
      "Skipping cell_line_efficacy.csv as it already exists, pass `overwrite=True` to overwrite\n",
      "Skipping library_effect.csv as it already exists, pass `overwrite=True` to overwrite\n",
      "Skipping t0_offset.csv as it already exists, pass `overwrite=True` to overwrite\n",
      "all files fetched, tranforming format\n",
      "transforming guide efficacy\n",
      "done\n"
     ]
    }
   ],
   "source": [
    "chronos.fetch_parameters()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create a model with the data we want to train as before, but with two important details:\n",
    "- we pass the argument `pretrained=True` when we initialize\n",
    "- the library batch names must match the DepMap library batch names, as that's what we're using for the pretrained model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normalizing readcounts\n",
      "\n",
      "\n",
      "Finding all unique guides and genes\n",
      "found 3474 unique guides and 883 unique genes in Achilles-Avana-2D\n",
      "found 4084 unique guides and 833 unique genes in Achilles-KY-2D\n",
      "found 7558 unique guides and 887 unique genes overall\n",
      "\n",
      "finding guide-gene mapping indices\n",
      "\n",
      "finding all unique sequenced replicates, cell lines, and pDNA batches\n",
      "found 92 unique sequences (excluding pDNA) and 44 unique cell lines in Achilles-Avana-2D\n",
      "found 63 unique sequences (excluding pDNA) and 23 unique cell lines in Achilles-KY-2D\n",
      "found 155 unique replicates and 58 unique cell lines overall\n",
      "\n",
      "finding replicate-cell line mappings indices\n",
      "\n",
      "finding replicate-pDNA mappings indices\n",
      "\n",
      "\n",
      "assigning float constants\n",
      "Estimating or aligning variances\n",
      "\tEstimating excess variance (alpha) for Achilles-Avana-2D\n",
      "\tEstimating excess variance (alpha) for Achilles-KY-2D\n",
      "Creating excess variance tensors\n",
      "\tCreated excess variance tensor for Achilles-Avana-2D with shape [92, 1]\n",
      "\tCreated excess variance tensor for Achilles-KY-2D with shape [63, 1]\n",
      "initializing graph\n",
      "\n",
      "building gene effect mask\n",
      "\n",
      "building doubling vectors\n",
      "made days vector of shape [92, 1] for Achilles-Avana-2D\n",
      "made days vector of shape [63, 1] for Achilles-KY-2D\n",
      "\n",
      "building late observed timepoints\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dempster/packages/chronos_public/chronos/model.py:377: UserWarning: Between 0 (batch=Index(['Achilles-Avana-2D_Avana-2', 'Achilles-Avana-2D_Avana-3',\n",
      "       'Achilles-Avana-2D_Avana-4'],\n",
      "      dtype='object')) and 0 (batch=Index(['Achilles-Avana-2D_Avana-2', 'Achilles-Avana-2D_Avana-3',\n",
      "       'Achilles-Avana-2D_Avana-4'],\n",
      "      dtype='object')) negative control sgRNAs were found to be systematically over- or under-represented in the screens and excluded.\n",
      "  warn(\"Between %i (batch=%r) and %i (batch=%r) negative control sgRNAs were found to be \\\n",
      "/Users/dempster/packages/chronos_public/chronos/model.py:377: UserWarning: Between 4 (batch=Index(['Achilles-KY-2D_KY-1'], dtype='object')) and 39 (batch=Index(['Achilles-KY-2D_KY-2'], dtype='object')) negative control sgRNAs were found to be systematically over- or under-represented in the screens and excluded.\n",
      "  warn(\"Between %i (batch=%r) and %i (batch=%r) negative control sgRNAs were found to be \\\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\tbuilt normalized timepoints for Achilles-Avana-2D with shape (92, 3474) (replicates X guides)\n",
      "\tbuilt normalized timepoints for Achilles-KY-2D with shape (63, 4084) (replicates X guides)\n",
      "\n",
      "building t0 reads\n",
      "\n",
      "\n",
      "Building variables\n",
      "\n",
      "building t0 reads estimate\n",
      "made t0 batch with shape (2, 4084) for Achilles-Avana-2D\n",
      "made t0 batch with shape (2, 4084) for Achilles-KY-2D\n",
      "building guide efficacy\n",
      "built guide efficacy: shape [1, 7559]\n",
      "building growth rate\n",
      "built growth rate: shape {'Achilles-Avana-2D': [58, 1], 'Achilles-KY-2D': [58, 1]}\n",
      "\n",
      "building cell line efficacy\n",
      "built cell line efficacy: shapes {'Achilles-Avana-2D': [58, 1], 'Achilles-KY-2D': [58, 1]}\n",
      "building screen delay\n",
      "built screen delay\n",
      "building gene effect\n",
      "built core gene effect: 58 cell lines by 887 genes\n",
      "\n",
      "\n",
      "Connecting graph nodes in model\n",
      "\n",
      "building effective days\n",
      "built effective days, shapes {'Achilles-Avana-2D': [92, 887], 'Achilles-KY-2D': [63, 887]}\n",
      "\n",
      "building gene effect growth graph nodes\n",
      "built gene effect growth graph nodes, shapes {'Achilles-Avana-2D': [58, 887], 'Achilles-KY-2D': [58, 887]}\n",
      "\n",
      "building combined efficacy\n",
      "built combined efficacy, shape {'Achilles-Avana-2D': [58, 3474], 'Achilles-KY-2D': [58, 4084]}\n",
      "built expanded combined efficacy, shapes {'Achilles-Avana-2D': [92, 3474], 'Achilles-KY-2D': [63, 4084]}\n",
      "\n",
      "building growth estimates of edited cells and overall estimates of fold change in guide abundance\n",
      "built growth and change\n",
      "\n",
      "building unnormalized estimates of final abundance\n",
      "built unnormalized abundance\n",
      "\n",
      "building normalized estimates of final abundance\n",
      "built normalized abundance\n",
      "initializing precost variables\n",
      "\n",
      "\n",
      "Building all costs\n",
      "\n",
      "assembling guide efficacy regularization\n",
      "building smoothed regularization\n",
      "\n",
      "building t0 reads regularization/cost\n",
      "\n",
      "building NB2 cost\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:57:47.976130: W tensorflow/c/c_api.cc:304] Operation '{name:'inferred_t0_3/base_Achilles-Avana-2D/Assign' id:11469 op device:{requested: '', assigned: ''} def:{{{node inferred_t0_3/base_Achilles-Avana-2D/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](inferred_t0_3/base_Achilles-Avana-2D, inferred_t0_3/base_Achilles-Avana-2D/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "building other regularizations\n",
      "\n",
      "Creating optimizer\n",
      "initializing rest of graph\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-26 13:57:49.364178: W tensorflow/c/c_api.cc:304] Operation '{name:'GE_3/library_effect/Achilles-KY-2D/Assign' id:11677 op device:{requested: '', assigned: ''} def:{{{node GE_3/library_effect/Achilles-KY-2D/Assign}} = AssignVariableOp[_has_manual_control_dependencies=true, dtype=DT_DOUBLE, validate_shape=false](GE_3/library_effect/Achilles-KY-2D, GE_3/library_effect/Achilles-KY-2D/Initializer/initial_value)}}' was changed by setting attribute after it was run by a session. This mutation will have no effect, and will trigger an error in the future. Either don't modify nodes after running them or create a new session.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimating initial screen efficacy and gene effect\n",
      "\t Achilles-Avana-2D\n",
      "\t Achilles-KY-2D\n",
      "\tverifying graph integrity\n",
      "verifying user inputs\n",
      "verifying variables\n",
      "verifying calculated terms\n",
      "\tAchilles-Avana-2D _gene_effect\n",
      "\tAchilles-Avana-2D _selected_efficacies\n",
      "\tAchilles-Avana-2D_predicted_readcounts_unscaled\n",
      "\tAchilles-Avana-2D _predicted_readcounts\n",
      "\tAchilles-Avana-2D _normalized_readcounts\n",
      "\tAchilles-Avana-2D _cost_presum\n",
      "sess run\n",
      "\tAchilles-Avana-2D _cost\n",
      "\tAchilles-Avana-2D _full_costs\n",
      "\tAchilles-KY-2D _gene_effect\n",
      "\tAchilles-KY-2D _selected_efficacies\n",
      "\tAchilles-KY-2D_predicted_readcounts_unscaled\n",
      "\tAchilles-KY-2D _predicted_readcounts\n",
      "\tAchilles-KY-2D _normalized_readcounts\n",
      "\tAchilles-KY-2D _cost_presum\n",
      "sess run\n",
      "\tAchilles-KY-2D _cost\n",
      "\tAchilles-KY-2D _full_costs\n",
      "waiting for user to load model\n"
     ]
    }
   ],
   "source": [
    "model2_pretrained = chronos.Chronos(\n",
    "    sequence_map={\"Achilles-Avana-2D\": sequence_map, 'Achilles-KY-2D': ky_sequence_map},\n",
    "    guide_gene_map={\"Achilles-Avana-2D\": guide_map, 'Achilles-KY-2D': ky_guide_map},\n",
    "    readcounts={\"Achilles-Avana-2D\": readcounts, 'Achilles-KY-2D': ky_readcounts},\n",
    "    negative_control_sgrnas={\"Achilles-Avana-2D\": negative_controls, \"Achilles-KY-2D\": ky_negative_controls},\n",
    "    pretrained=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we import the DepMap data from the directory into the model, and train:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chronos model loaded\n",
      "ready to train\n"
     ]
    }
   ],
   "source": [
    "model2_pretrained.import_model(\"./Data/DepMapParameters/\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NB2 cost 0.27116384868792104\n",
      "Full cost 0.46402864438737385\n",
      "relative_growth_rate\n",
      "\tAchilles-Avana-2D max 1.039, min 0.95394\n",
      "\tAchilles-KY-2D max 1.047, min 0.96608\n",
      "mean guide efficacy 0.8814048828522121\n",
      "t0_offset SD: [('Achilles-Avana-2D', 0.1755915504675398), ('Achilles-KY-2D', 0.14386042441508556)]\n",
      "\n",
      "gene mean -0.44947516580140634\n",
      "SD of gene means 0.6593218864040892\n",
      "Mean of gene SDs 0.3420810837838141\n",
      "\n",
      "\n",
      "\n",
      "51 epochs trained, time taken 0:00:01, projected remaining 0:00:07\n",
      "NB2 cost 0.24021191845091294\n",
      "Full cost 0.34959317695282827\n",
      "relative_growth_rate\n",
      "\tAchilles-Avana-2D max 1.476, min 0.55084\n",
      "\tAchilles-KY-2D max 1.363, min 0.63594\n",
      "mean guide efficacy 0.8814048828522121\n",
      "t0_offset SD: [('Achilles-Avana-2D', 0.1755915504675398), ('Achilles-KY-2D', 0.14386042441508556)]\n",
      "\n",
      "gene mean -0.3992848841580951\n",
      "SD of gene means 0.6932161219723673\n",
      "Mean of gene SDs 0.28846565769376614\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model2_pretrained.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3.11",
   "language": "python",
   "name": "python3.11"
  },
  "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.11.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "183.991px",
    "width": "287.997px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "272px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}