cell-cycle/Results.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results of the Cytometry experiment\n",
    "\n",
    "### Summary of the experiment\n",
    "\n",
    "1. Data gathering from [BBBC022](https://data.broadinstitute.org/bbbc/BBBC022/) \n",
    "\n",
    "2. Compound selection and data conditioning (`Compound_Images_List.ipynb`)\n",
    "\n",
    "3. Segmentation and DNA content analysis (for UNET)\n",
    "    * Segmentation (`batch_prediction.py`)\n",
    "    * DNA content analysis\n",
    "\n",
    "4. Segmentation and DNA content analysis (for CellProfiler)\n",
    "    * Segmentation\n",
    "    * DNA content analysis\n",
    "\n",
    "5. Results exploitation (`Results.ipynb`)\n",
    "\n",
    "\n",
    "### 1. Data gathering\n",
    "\n",
    "The dataset is coming from the [Broad Bioimage Benchmark Collection](https://data.broadinstitute.org/bbbc/index.html) and its name is [BBBC022](https://data.broadinstitute.org/bbbc/BBBC022/). The channel we are interested in is the Hoechst 33342 channel, so we only need to download the `.zip` files ending with `w1`. All the files must be decompressed and stored in a `/data` folder.\n",
    "\n",
    "\n",
    "### 2. Compound selection and data conditioning\n",
    "\n",
    "The list of the available compounds is printed by the first cells of `Compound_Images_List.ipynb`. In the section \"Compound selection\" we can choose which compound to analyse.\n",
    "\n",
    "Once the compound is selected, we can run the entire notebook `Compound_Images_List.ipynb` for the desired compound. It should output a `.csv` file in the current folder with the list of the images corresponding to the compound. It should also create the structure of the directories in order to use the `batch_prediction.py` code.\n",
    "\n",
    "It creates a `data/COMPOUND` directory with three subdirectories : \n",
    "* `/COMPOUND/images` : folder of the images corresponding to the compound (filled by `Compound_Images_List.ipynb`)\n",
    "* `/COMPOUND/masks` : folder of the masks found with UNET (computed by `batch_prediction.py`)\n",
    "* `/COMPOUND/cp` : folder of the intermediary output of the DNA content analysis with CellProfiler\n",
    "\n",
    "It also creates a `/COMPOUND` folder in another location where the results of the CellProfiler segmentation will be stored.\n",
    "\n",
    "\n",
    "### 3. Segmentation and DNA content analysis (for UNET)\n",
    "\n",
    "#### 3. a) Segmentation\n",
    "\n",
    "Use the `batch_prediction.py` file to process the images. It should be called by :\n",
    "\n",
    "```python batch_prediction.py experiment_name image_list.csv input_dir output_dir```\n",
    "\n",
    "The recommended structure is :\n",
    "\n",
    "```python3 batch_prediction.py experiment1 ./../../DSB\\ paper/COMPOUND_image_list.csv ./../../UNET/COMPOUND/images/ ./../../UNET/COMPOUND/masks/```\n",
    "\n",
    "where `COMPOUND` is the selected compound for the compound analysis.\n",
    "\n",
    "\n",
    "*Warnings:* \n",
    "* This command should be called in the folder where `batch_prediction.py` is\n",
    "* python3 can be used instead of python\n",
    "* The experiment name can be chosen by the user\n",
    "* The list of images is the output of the compound analysis\n",
    "* The `model.hdf5` file should be located in the `experiment/experiment_name` folder\n",
    "\n",
    "This step takes care of the segmentation using UNET. It outputs the mask corresponding of each image in `COMPOUND_image_list.csv `.\n",
    "\n",
    "\n",
    "#### 3. b) DNA content analysis\n",
    "\n",
    "Use a pipeline in CellProfiler to compute the integrated intensity of the images, masked by the segmentation masks from the prediction. The pipeline consists of a `MeasureObjectIntensity` layer and some other layers to save the results.\n",
    "\n",
    "For the several compounds (Gabazine, CYCLOPIAZONIC ACID, ESTRADIOL, DMSO, Proglumide, PIZOTIFEN, CATECHIN, PENTABENZOATE, TRIADIMEFON, 3-HYDROXYCOUMARIN, Etoposide, DACTINOMYCIN, Colchicine, Blebbistatin) the projects are saved into their respective folders (`/data/COMPOUND`). The pipeline for another compound should have exactly the same structure in order to use this notebook. \n",
    "\n",
    "\n",
    "### 4. Segmentation and DNA content analysis (for CellProfiler)\n",
    "\n",
    "The segmentation and the DNA content analysis for CellProfiler are done using the same pipeline. The pipeline consists of one step of segmentation (`IdentifyPrimaryObjects` layer) and one step of measurement (`MeasureObjectIntensity`) to compute the DNA content of the objects.\n",
    "\n",
    "\n",
    "### 5. Results exploitation\n",
    "\n",
    "DNA content analysis is the purpose of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "\n",
    "import math\n",
    "\n",
    "import skimage\n",
    "from skimage import io\n",
    "\n",
    "from os import listdir\n",
    "from os.path import isfile, join"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results gathering\n",
    "\n",
    "\n",
    "The following functions load the `.csv` file corresponding to the DNA content analysis for UNET and CellProfiler. If the above structure is respected, they should work perfectly, otherwise the paths may need to be adapted.\n",
    "\n",
    "They create for each segmentation method (UNET and CP) a dataframe containing the information of the analysis :\n",
    "\n",
    "| Index | ImageNumber | ObjectNumber | Integrated_Intensity | Plate |\n",
    "|:-----:|-------------|--------------|----------------------|-------|\n",
    "| **0** | 1           | 1            | 2.245075             | 20589 |\n",
    "| **1** | 1           | 2            | 1.579263             | 20589 |\n",
    "| ...   | ...         | ...          | ...                  | ...   |\n",
    "\n",
    "\n",
    "Each *ImageNumber* corresponds to an image, each *ObjectNumber* corresponds to one found nucleus. The *Integrated_Intensity* field shows the result of the DNA content analysis, and we keep the information of the plate.\n",
    "\n",
    "For each compound analysis, we have 4 dataframes of the above type :\n",
    "\n",
    "|  Compound\\Method | UNET         | CP         |\n",
    "|:----------------:|--------------|------------|\n",
    "| **Compound X :** | df_unet_cmpd | df_cp_cmpd |\n",
    "|       **DMSO :** | df_unet_dmso | df_cp_dmso |\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_dataframe_unet(data_path, compound_number):\n",
    "    compounds = ['Gabazine', 'CYCLOPIAZONIC ACID', \n",
    "                 'ESTRADIOL', 'DMSO', \n",
    "                 'Proglumide', 'PIZOTIFEN', \n",
    "                 'CATECHIN PENTABENZOATE', 'TRIADIMEFON',\n",
    "                 '3-HYDROXYCOUMARIN', 'Etoposide',\n",
    "                 'DACTINOMYCIN', 'Colchicine',\n",
    "                 'Blebbistatin']\n",
    "\n",
    "    compounds2 = ['Gabazine', 'CYCLOPIAZONIC_ACID', \n",
    "                  'ESTRADIOL', 'DMSO', \n",
    "                  'Proglumide', 'PIZOTIFEN', \n",
    "                  'CATECHIN_PENTABENZOATE', 'TRIADIMEFON',\n",
    "                  '3-HYDROXYCOUMARIN', 'Etoposide',\n",
    "                  'DACTINOMYCIN', 'Colchicine',\n",
    "                  'Blebbistatin']\n",
    "\n",
    "    compound = compounds[compound_number]\n",
    "    compound2 = compounds2[compound_number]\n",
    "\n",
    "    mypath = data_path + compound\n",
    "    onlyfiles = [f for f in listdir(mypath) if isfile(join(mypath, f)) and f not in ['UNET_MASKS.csv', 'UNET_Image.csv']]\n",
    "    masks_end = 'MASKS.csv'\n",
    "    image_end = 'Image.csv'\n",
    "    masks_file = mypath + '/'\n",
    "    image_file = mypath + '/'\n",
    "\n",
    "    for f in onlyfiles:\n",
    "        if (len(f) < len(masks_end) or len(f) < len(image_end)):\n",
    "            continue\n",
    "\n",
    "        if (f[len(f)-len(masks_end):] == masks_end):\n",
    "            masks_file += f\n",
    "\n",
    "        if (f[len(f)-len(image_end):] == image_end):\n",
    "            image_file += f\n",
    "\n",
    "    df_masks = pd.read_csv(masks_file)\n",
    "    df_image = pd.read_csv(image_file)\n",
    "\n",
    "    df1 = df_masks[['ImageNumber', 'ObjectNumber','Intensity_IntegratedIntensity_'+compound2+'_IMAGES']]\n",
    "    df2 = df_image[['FileName_'+compound2+'_IMAGES']]\n",
    "\n",
    "    ser = []\n",
    "    for i in range(0, len(df1)):\n",
    "        ser.append(df2['FileName_'+compound2+'_IMAGES'][df1.ImageNumber.iloc[i]-1][18:23])\n",
    "\n",
    "    ser = np.array(ser)\n",
    "    df1 = df1.assign(Plate=pd.Series(ser).values)\n",
    "\n",
    "    df1 = df1.rename(index=str, columns={'Intensity_IntegratedIntensity_'+compound2+'_IMAGES': \"Integrated_Intensity\"})\n",
    "\n",
    "    return df1, compound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_dataframe_unet2(data_path, compound_number):\n",
    "    compounds = ['Gabazine', 'CYCLOPIAZONIC ACID', \n",
    "                 'ESTRADIOL', 'DMSO', \n",
    "                 'Proglumide', 'PIZOTIFEN', \n",
    "                 'CATECHIN PENTABENZOATE', 'TRIADIMEFON',\n",
    "                 '3-HYDROXYCOUMARIN', 'Etoposide',\n",
    "                 'DACTINOMYCIN', 'Colchicine',\n",
    "                 'Blebbistatin']\n",
    "\n",
    "    compound = compounds[compound_number]\n",
    "\n",
    "    mypath = data_path + compound\n",
    "    masks_end = 'MASKS.csv'\n",
    "    image_end = 'Image.csv'\n",
    "    masks_file = mypath + '/UNET_DilateObjects.csv'\n",
    "    image_file = mypath + '/UNET_Image.csv'\n",
    "    \n",
    "    df_masks = pd.read_csv(masks_file)\n",
    "    df_image = pd.read_csv(image_file)\n",
    "\n",
    "    df1 = df_masks[['ImageNumber', 'ObjectNumber','Intensity_IntegratedIntensity_IMAGES']]\n",
    "    df2 = df_image[['FileName_IMAGES']]\n",
    "\n",
    "    ser = []\n",
    "    for i in range(0, len(df1)):\n",
    "        ser.append(df2['FileName_IMAGES'][df1.ImageNumber.iloc[i]-1][18:23])\n",
    "\n",
    "    ser = np.array(ser)\n",
    "    df1 = df1.assign(Plate=pd.Series(ser).values)\n",
    "\n",
    "    df1 = df1.rename(index=str, columns={'Intensity_IntegratedIntensity_IMAGES': \"Integrated_Intensity\"})\n",
    "\n",
    "    return df1, compound"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def construct_dataframe_cp(data_path, compound_number):\n",
    "    compounds = ['Gabazine', 'CYCLOPIAZONIC ACID', \n",
    "                 'ESTRADIOL', 'DMSO', \n",
    "                 'Proglumide', 'PIZOTIFEN', \n",
    "                 'CATECHIN PENTABENZOATE', 'TRIADIMEFON',\n",
    "                 '3-HYDROXYCOUMARIN', 'Etoposide',\n",
    "                 'DACTINOMYCIN', 'Colchicine',\n",
    "                 'Blebbistatin']\n",
    "    compound = compounds[compound_number]\n",
    "    compounds2 = ['Gabazine', 'CYCLOPIAZONIC_ACID', \n",
    "                  'ESTRADIOL', 'DMSO', \n",
    "                  'Proglumide', 'PIZOTIFEN', \n",
    "                  'CATECHIN_PENTABENZOATE', 'TRIADIMEFON',\n",
    "                  '3-HYDROXYCOUMARIN', 'Etoposide',\n",
    "                  'DACTINOMYCIN', 'Colchicine',\n",
    "                  'Blebbistatin']\n",
    "    compound2 = compounds2[compound_number]\n",
    "\n",
    "    compound_path = data_path + compound + '/'\n",
    "    \n",
    "    df1 = pd.read_csv(compound_path + 'Nuclei.csv')\n",
    "    df2 = pd.read_csv(compound_path + 'Image.csv')\n",
    "\n",
    "    df1 = df1[['ImageNumber', 'ObjectNumber','Intensity_IntegratedIntensity_DNA']]\n",
    "    df2 = df2[['FileName_DNA']]\n",
    "    \n",
    "    ser = []\n",
    "    for i in range(0, len(df1)):\n",
    "        ser.append(df2['FileName_DNA'][df1.ImageNumber.iloc[i]-1][18:23])\n",
    "\n",
    "    ser = np.array(ser)\n",
    "    df1 = df1.assign(Plate=pd.Series(ser).values)\n",
    "\n",
    "    df1 = df1.rename(index=str, columns={'Intensity_IntegratedIntensity_DNA': \"Integrated_Intensity\"})\n",
    "    \n",
    "    return df1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Plotting \n",
    "\n",
    "For each of the 2 methods and each compound, one can print an image, its corresponding mask, and the image seen through the mask thanks to the 2 functions below :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_image_and_mask_unet(data_path, compound, i):\n",
    "    mypath = data_path + compound\n",
    "    \n",
    "    image_path = mypath + '/images/'\n",
    "    masks_path = mypath + '/masks/'\n",
    "    \n",
    "    onlyfiles = [f for f in listdir(image_path) if isfile(join(image_path, f))]\n",
    "    im_end = '.tif'\n",
    "    files = [f for f in onlyfiles if f[-4:]==im_end]\n",
    "    \n",
    "    file = files[i]\n",
    "#     print(file)\n",
    "    \n",
    "    ipath = image_path + file\n",
    "    mpath = masks_path + file\n",
    "    \n",
    "    im = skimage.io.imread(ipath)\n",
    "    ma = skimage.io.imread(mpath)\n",
    "\n",
    "    plt.figure(figsize=(15,10))\n",
    "    plt.subplot(1,2,1)\n",
    "    plt.imshow(im)\n",
    "    plt.subplot(1,2,2)\n",
    "#     plt.imshow(ma>0)\n",
    "    plt.imshow(ma)\n",
    "    \n",
    "    plt.figure(figsize=(7,6))\n",
    "    plt.imshow(im*(ma>0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_image_and_mask_cp(data_path_im, data_path_cp, compound, i):\n",
    "    mypath = data_path_im + compound\n",
    "    \n",
    "    image_path = mypath + '/images/'\n",
    "    \n",
    "    onlyfiles = [f for f in listdir(image_path) if isfile(join(image_path, f))]\n",
    "    im_end = '.tif'\n",
    "    files = [f for f in onlyfiles if f[-4:]==im_end]\n",
    "    \n",
    "    file = files[i]\n",
    "#     print(file)\n",
    "    ipath = image_path + file\n",
    "    \n",
    "    \n",
    "    mpath = data_path_cp + compound +'/' + file[:-4] + '_Overlay.tiff'\n",
    "#     print(mpath)\n",
    "    \n",
    "    im = skimage.io.imread(ipath)\n",
    "    ma = skimage.io.imread(mpath)\n",
    "    plt.figure(figsize=(15,10))\n",
    "    plt.subplot(1,2,1)\n",
    "    plt.imshow(im)\n",
    "    plt.subplot(1,2,2)\n",
    "    plt.imshow(ma)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Normalization\n",
    "\n",
    "Each plate corresponds to different experimental conditions. Hence, we cannot compare the intensity of a cell in plate A with the intensity of a cell in plate B without normalizing them.\n",
    "\n",
    "To normalize the intensities over the different plates, we use the log-transform of the histograms. \n",
    "\n",
    "\n",
    "### Splitting and log-intensities\n",
    "\n",
    "We first split the dataframe into a dictionary with each plate as entry (`split_per_plates`) and then, we compute the log10 of each intensity value (`compute_log_intensities`), for each plate.\n",
    "\n",
    "This gives us a dictionary of dataframes :\n",
    "\n",
    "{'20589':\n",
    "\n",
    "| Index | ImageNumber | ObjectNumber | Integrated_Intensity | Plate | Log_Integrated_Intensity |\n",
    "|:-----:|-------------|--------------|----------------------|-------|--------------------------|\n",
    "| **0** | 1           | 1            | 2.245075             | 20589 | 0.351231                 |\n",
    "| **1** | 1           | 2            | 1.579263             | 20589 | 0.198454                 |\n",
    "| ...   | ...         | ...          | ...                  | ...   | ...                      |\n",
    "\n",
    "'20590':\n",
    "\n",
    "| Index | ImageNumber | ObjectNumber | Integrated_Intensity | Plate | Log_Integrated_Intensity |\n",
    "|:-----:|-------------|--------------|----------------------|-------|--------------------------|\n",
    "| **0** | 10          | 1            | 3.067079             | 20590 | 0.486725                 |\n",
    "| **1** | 10          | 2            | 2.338338             | 20590 | 0.368907                 |\n",
    "| ...   | ...         | ...          | ...                  | ...   | ...                      |\n",
    "\n",
    "...}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_per_plates(df):\n",
    "    plates = df.Plate.unique()\n",
    "    dfs = {}\n",
    "\n",
    "    for p in plates:\n",
    "        dfs[p] = df[df.Plate == p]\n",
    "    \n",
    "    return dfs\n",
    "\n",
    "def compute_log_intensities(dfs):\n",
    "    for p in dfs:\n",
    "        LII = np.log10(dfs[p].Integrated_Intensity.values)\n",
    "        dfs[p] = dfs[p].assign(Log_Integrated_Intensity = pd.Series(LII).values)\n",
    "    return dfs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization\n",
    "\n",
    "To normalize the instensities over the plates, we use the log intensities. If we look at the typical DNA content distribution, it should contain the global maximum at position 2N and a local maximum at position 4N. This is exactly what we obtain for each plate, but the position 2N is not the same for all of the plates (thus the position 4N is not neither). By shifting the histogram in the log space, we can easily make the two modes match in the original space. That is why we use the log transform.\n",
    "\n",
    "Hence, in order to normalize, we compute the histogram of the log intensities. The goal is to find the main peak for each plates and to align these peaks. To find the main peak for each plate, we first smoothen the histogram to avoid the influence of the noise, and we take the maximum.\n",
    "\n",
    "Once we computed the position of the maximum for each plate, we need to align them so the main peak of each histogram is at the same position. To do so, we take one of the positions and we compute the shifts needed to reach this position from the others. We can now apply those shifts to the log intensities (each shift corresponds to a plate).\n",
    "\n",
    "Since the compound could have an influence on the cell cycle, we take the DMSO to compute the shifts. Since same plate means same experimental conditions, the shift computed to normalize the DMSO on a plate can be used to normalize the compound on the same plate.\n",
    "\n",
    "\n",
    "Finally, at the end of the function, we remove the plates that are not in the compound dictionary from the dictionary of the DMSO to keep only the plates of interest.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalization(dfs_cmpd, dfs_dmso, plot_histograms=False, return_shifts=False):\n",
    "    dfs = dfs_dmso.copy()\n",
    "    dfs2 = dfs_cmpd.copy()\n",
    "    \n",
    "    nbins = 200\n",
    "    hrange = (-1.0, 2.5)\n",
    "    step = (hrange[1]-hrange[0])/nbins\n",
    "\n",
    "    n_mask = 7\n",
    "    mask = (1/n_mask)*np.ones(n_mask)\n",
    "\n",
    "    maxima = {}\n",
    "\n",
    "    for plate in dfs:\n",
    "        histogram = np.histogram(\n",
    "            a=dfs[plate].Log_Integrated_Intensity,\n",
    "            bins = nbins,\n",
    "            range = hrange\n",
    "        )\n",
    "        \n",
    "        sig = np.convolve(histogram[0], mask)\n",
    "        maxima[plate] = np.argmax(sig)\n",
    "\n",
    "    maxis = np.fromiter(maxima.values(), dtype=float)\n",
    "    max_pos = np.max(maxis)\n",
    "\n",
    "    shifts = {}\n",
    "    for plate in maxima:\n",
    "        shifts[plate] = (max_pos-maxima[plate])*step\n",
    "\n",
    "    for plate in dfs2:\n",
    "        vals = dfs2[plate].Log_Integrated_Intensity.values\n",
    "        dfs2[plate] = dfs2[plate].assign(normalized_logInt = pd.Series(vals+shifts[plate]).values)\n",
    "        dfs2[plate] = dfs2[plate].assign(Normalized_Intensity = pd.Series(np.power(10, dfs2[plate].normalized_logInt)).values)\n",
    "        dfs_cmpd[plate] = dfs2[plate].assign(Normalized_Intensity = pd.Series(dfs2[plate].Normalized_Intensity).values)\n",
    "        \n",
    "        vals = dfs[plate].Log_Integrated_Intensity.values\n",
    "        dfs[plate] = dfs[plate].assign(normalized_logInt = pd.Series(vals+shifts[plate]).values)\n",
    "        dfs[plate] = dfs[plate].assign(Normalized_Intensity = pd.Series(np.power(10, dfs[plate].normalized_logInt)).values)\n",
    "        dfs_dmso[plate] = dfs[plate].assign(Normalized_Intensity = pd.Series(dfs[plate].Normalized_Intensity).values)\n",
    "        \n",
    "            \n",
    "    if plot_histograms:\n",
    "        for plate in dfs_cmpd:\n",
    "            plt.figure(figsize=(20,8))\n",
    "            \n",
    "#             print(\"Plate {} : {}\".format(plate,len(dfs[plate])))\n",
    "            plt.subplot(1, 2, 1)\n",
    "            plt.title(plate + ' DMSO')\n",
    "            p1 = plt.hist(dfs_dmso[plate].Integrated_Intensity, nbins, (0, 30), label='Before normalization')\n",
    "            p2 = plt.hist(dfs_dmso[plate].Normalized_Intensity, nbins, (0, 30), label='After normalization')\n",
    "            plt.xlabel('Logarithm intensities')\n",
    "            plt.ylabel('number of cells')\n",
    "            plt.legend(frameon = False)\n",
    "            \n",
    "            plt.subplot(1, 2, 2)\n",
    "            plt.title(plate + ' COMPOUND')\n",
    "            plt.hist(dfs_cmpd[plate].Integrated_Intensity, nbins, (0, 30), label='Before normalization')\n",
    "            plt.hist(dfs_cmpd[plate].Normalized_Intensity, nbins, (0, 30), label='After normalization')\n",
    "            plt.xlabel('Intensities')\n",
    "            plt.ylabel('number of cells')\n",
    "            plt.legend(frameon = False)\n",
    "    \n",
    "    dfs_dmso2 = dfs_dmso.copy()\n",
    "    for plate in dfs_dmso2:\n",
    "        if plate not in dfs_cmpd.keys():\n",
    "            del dfs_dmso[plate]\n",
    "            \n",
    "    if return_shifts:\n",
    "        return dfs_cmpd, dfs_dmso, shifts\n",
    "    \n",
    "    return dfs_cmpd, dfs_dmso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regrouping the plates\n",
    "\n",
    "Now that all the plates are normalized, we need to merge them together in order to compute the Z' factor for all the images. That is the purpose of the function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regroup_plates(df, dfs):\n",
    "    df2 = df.copy()\n",
    "    pls = df2.Plate.unique()\n",
    "    pls2 = dfs.keys()\n",
    "    \n",
    "    plates = list(set(list(pls))-set(list(pls2)))\n",
    "    \n",
    "    for p in plates:\n",
    "        df2 = df2[df2.Plate != p]\n",
    "\n",
    "    intensities = []\n",
    "\n",
    "    for plate in pls2:\n",
    "        intensities = intensities + list(dfs[plate].Normalized_Intensity.values)\n",
    "    \n",
    "    \n",
    "    df2 = df2.assign(Normalized_Intensity=pd.Series(np.array(intensities)).values)\n",
    "    \n",
    "    return df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regroup_plates2(dfs):\n",
    "    plates = list(dfs.keys())\n",
    "    df = pd.DataFrame(columns=dfs[plates[0]].columns)\n",
    "\n",
    "    for plate in plates:\n",
    "        df = df.append(dfs[plate])\n",
    "    \n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Z factor calculation\n",
    "\n",
    "\n",
    "We compute the Z factor to determine if we can see a significant difference between the compound experiment results and the DMSO experiment results. \n",
    "\n",
    "\n",
    "### Percentages computation\n",
    "\n",
    "Thanks to the function below, we compute for each image the percentage of cell with intensity that is above a given threshold. By choosing the threshold between the 2 modes of the DNA distribution, we compute, for each image, the percentage of cells with intensity close to 4N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_percentages(df, threshold, display_hist=False, cmpd_name=''):\n",
    "    images = df.ImageNumber.unique()\n",
    "    percentages = []\n",
    "    for im in images:\n",
    "        vals = df[df.ImageNumber == im].Normalized_Intensity.values\n",
    "        percentages.append(np.sum(vals>threshold)/vals.shape[0])\n",
    "        \n",
    "    if display_hist:\n",
    "        image_number = 1\n",
    "        vals = np.array(percentages)\n",
    "    \n",
    "        x = np.array([threshold, threshold])\n",
    "        y = np.array([0, 100])\n",
    "        plt.figure(figsize=(10,8))\n",
    "        plt.title(compound + ' percentages')\n",
    "        plt.hist(vals, 100, (0.0, 1.0), label='Values', normed=True)\n",
    "        \n",
    "#         plt.plot(x, y, 'r', label='Threshold')\n",
    "        plt.xlabel('Proportion of intensity values above threshold')\n",
    "        plt.ylabel('Number of images')\n",
    "        plt.legend(frameon = False)\n",
    "        \n",
    "    return percentages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Threshold computation\n",
    "\n",
    "To compute the threshold, we first find the position of the peak at 2N and 4N, and we take the mean between the two positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_threshold(df_DMSO, display=False):\n",
    "    nbins = 200\n",
    "    hrange = (-1.0, 30)\n",
    "    step = (hrange[1]-hrange[0])/nbins\n",
    "\n",
    "    n_mask = 3\n",
    "    mask = (1/n_mask)*np.ones(n_mask)\n",
    "\n",
    "    histogram2 = np.histogram(\n",
    "        a=df_DMSO.Normalized_Intensity,\n",
    "        bins = nbins,\n",
    "        range = hrange\n",
    "    )\n",
    "    sig = np.convolve(histogram2[0], mask)\n",
    "    maxima = np.argmax(sig)*step\n",
    "    \n",
    "    local_maxima = np.r_[True, sig[1:] > sig[:-1]] & np.r_[sig[:-1] > sig[1:], True]\n",
    "    pos_max = np.squeeze(np.argwhere(local_maxima==True))\n",
    "    val_max = sig[local_maxima==True]\n",
    "    \n",
    "    argmax1 = np.argmax(val_max)\n",
    "    pmax1 = pos_max[argmax1]\n",
    "    max1 = val_max[argmax1]\n",
    "    \n",
    "    pos_max = pos_max[argmax1+1:]\n",
    "    val_max = val_max[argmax1+1:]\n",
    "    \n",
    "    argmax2 = np.argmax(val_max)\n",
    "    pmax2 = pos_max[argmax2]\n",
    "    \n",
    "    if display:\n",
    "        height = np.array([0, 1100])\n",
    "        plt.figure()\n",
    "        plt.plot(np.array(range(-1, nbins+1))*step, sig)\n",
    "        plt.plot(np.array([pmax1, pmax1])*step, height, 'r')\n",
    "        plt.plot(np.array([pmax2, pmax2])*step, height, 'r')\n",
    "    \n",
    "    thres = (pmax1+pmax2)/2*step\n",
    "    if display:\n",
    "        print(\"Threshold = {}\".format(thres))\n",
    "    \n",
    "    return thres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Z factor computation\n",
    "\n",
    "\n",
    "The formula of the Z factor is :\n",
    "\n",
    "$$Z' = 1-{3(\\sigma _{p}+\\sigma _{n}) \\over |\\mu _{p}-\\mu _{n}|}$$\n",
    "\n",
    "where p, n denote the positive and negative control, respectively compound and DMSO. $\\mu$ and $\\sigma$ are the mean and the standard deviation of the percentages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_Z_factor(percentages_pos, percentages_neg):\n",
    "    p = np.array(percentages_pos)\n",
    "    n = np.array(percentages_neg)\n",
    "    \n",
    "    mean_p = np.mean(p)\n",
    "    std_p  = np.std(p)\n",
    "    mean_n = np.mean(n)\n",
    "    std_n  = np.std(n)\n",
    "    \n",
    "    return 1 - 3*(std_p+std_n)/abs(mean_p-mean_n)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compound : Gabazine\n",
      "Z' for UNET = -28.066965940664478\n",
      "Z' for CellProfiler = -45.06801417449648\n",
      "\n",
      "Compound : CYCLOPIAZONIC ACID\n",
      "Z' for UNET = -29.996334244753626\n",
      "Z' for CellProfiler = -55.24586616428868\n",
      "\n",
      "Compound : ESTRADIOL\n",
      "Z' for UNET = -3.2878556884273085\n",
      "Z' for CellProfiler = -5.741775063390136\n",
      "\n",
      "Compound : Proglumide\n",
      "Z' for UNET = -9.799619561853978\n",
      "Z' for CellProfiler = -17.22519559743383\n",
      "\n",
      "Compound : PIZOTIFEN\n",
      "Z' for UNET = -39.49409987195662\n",
      "Z' for CellProfiler = -56.720938485163025\n",
      "\n",
      "Compound : CATECHIN PENTABENZOATE\n",
      "Z' for UNET = -59.539427374904186\n",
      "Z' for CellProfiler = -18.167502634304967\n",
      "\n",
      "Compound : TRIADIMEFON\n",
      "Z' for UNET = -73.22188573510053\n",
      "Z' for CellProfiler = -23.432665810567336\n",
      "\n",
      "Compound : 3-HYDROXYCOUMARIN\n",
      "Z' for UNET = -16.548270650652213\n",
      "Z' for CellProfiler = -51.00284333300692\n",
      "\n",
      "Compound : Etoposide\n",
      "Z' for UNET = -1.850535951808781\n",
      "Z' for CellProfiler = -1.7130831419201842\n",
      "\n",
      "Compound : DACTINOMYCIN\n",
      "Z' for UNET = -23.764863521752513\n",
      "Z' for CellProfiler = -301.046409211808\n",
      "\n",
      "Compound : Colchicine\n",
      "Z' for UNET = -13.48708224403214\n",
      "Z' for CellProfiler = -0.3252778152575291\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for cmp_number in range(0, 12):\n",
    "    if cmp_number == 3:\n",
    "        continue\n",
    "    \n",
    "    data_path = './../UNET/'\n",
    "\n",
    "    df_unet_cmpd, compound = construct_dataframe_unet2(data_path, cmp_number)\n",
    "    df_unet_dmso, dmso = construct_dataframe_unet2(data_path, 3)\n",
    "\n",
    "\n",
    "    df_plates_unet_cmpd, df_plates_unet_dmso = normalization(\n",
    "        compute_log_intensities(split_per_plates(df_unet_cmpd)),\n",
    "        compute_log_intensities(split_per_plates(df_unet_dmso)) , \n",
    "        plot_histograms=False)\n",
    "\n",
    "    df_unet_cmpd = regroup_plates2(df_plates_unet_cmpd)\n",
    "    df_unet_dmso = regroup_plates2(df_plates_unet_dmso)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    data_path_cp = './../CP/'\n",
    "\n",
    "    df_cp_cmpd = construct_dataframe_cp(data_path_cp, cmp_number)\n",
    "    df_cp_dmso = construct_dataframe_cp(data_path_cp, 3)\n",
    "\n",
    "    df_plates_cp_cmpd, df_plates_cp_dmso = normalization(\n",
    "        compute_log_intensities(split_per_plates(df_cp_cmpd)),\n",
    "        compute_log_intensities(split_per_plates(df_cp_dmso)),\n",
    "        plot_histograms=False)\n",
    "\n",
    "    df_cp_cmpd = regroup_plates2(df_plates_cp_cmpd)\n",
    "    df_cp_dmso = regroup_plates2(df_plates_cp_dmso)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    threshold_unet = compute_threshold(df_unet_dmso, display=False)\n",
    "    threshold_cp = compute_threshold(df_cp_dmso, display=False)\n",
    "\n",
    "    per_unet_cmpd = compute_percentages(df_unet_cmpd, threshold_unet)\n",
    "    per_unet_dmso = compute_percentages(df_unet_dmso, threshold_unet)\n",
    "\n",
    "    per_cp_cmpd = compute_percentages(df_cp_cmpd, threshold_cp)\n",
    "    per_cp_dmso = compute_percentages(df_cp_dmso, threshold_cp)\n",
    "\n",
    "    ZprimeUNET = compute_Z_factor(per_unet_cmpd, per_unet_dmso)\n",
    "    ZprimeCP = compute_Z_factor(per_cp_cmpd, per_cp_dmso)\n",
    "\n",
    "\n",
    "    print(\"Compound : {}\".format(compound))\n",
    "    print(\"Z' for UNET = {}\".format(ZprimeUNET))\n",
    "    print(\"Z' for CellProfiler = {}\".format(ZprimeCP))\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}