Image_Analysis_Final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "331126e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage import io, filters\n",
    "import aicsimageio\n",
    "from aicsimageio import AICSImage\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from matplotlib.patches import PathPatch\n",
    "from matplotlib.path import Path\n",
    "import matplotlib.patches as patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5614a50",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Function to plot protein expression\n",
    "#File is the AICS file containing the data\n",
    "#sig1 is the sigma value for the Gaussian filter for dapi\n",
    "#sig2 is the sigma value for the Gaussian filter for TdTomato\n",
    "#cut1 is the expression cutoff for dapi\n",
    "#cut2 is the expression cutoff for TdTomato\n",
    "#r1 and r2 are 2 element lists that define the rectangular ROI\n",
    "#One or both can also be zero if you do not want to define boundaries for that part\n",
    "#back1 and back2 are used to define a suitable background region containing no dapi signal\n",
    "\n",
    "def plot(file, sig1, sig2, cut1, cut2, r1, r2, back1, back2, plot = True):\n",
    "    #Read in the image and extract the channels\n",
    "    img = AICSImage(file)\n",
    "    goi = img.get_image_data(\"ZYX\", C=0, S=0, T=0)[0]\n",
    "    tdt = img.get_image_data(\"ZYX\", C=1, S=0, T=0)[0]\n",
    "    nes = img.get_image_data(\"ZYX\", C=2, S=0, T=0)[0]\n",
    "    dapi = img.get_image_data(\"ZYX\", C=3, S=0, T=0)[0]\n",
    "    \n",
    "    #If we want to plot, make a plot of the image including the ROI\n",
    "    if plot:\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.imshow(tdt)\n",
    "        rect = patches.Rectangle((r2[0], r1[0]), r2[1] - r2[0], r1[1] - r1[0], linewidth=1, edgecolor='r', facecolor='none')\n",
    "        ax.add_patch(rect)\n",
    "        plt.show()\n",
    "    \n",
    "    #plt.imshow(tdt)\n",
    "    #plt.show()\n",
    "    #plt.imshow(goi)\n",
    "    #plt.show()\n",
    "    #plt.imshow(dapi)\n",
    "    #plt.show()\n",
    "    \n",
    "    #Apply Gaussian filter\n",
    "    dapiFilt = filters.gaussian(dapi, sigma=sig1)\n",
    "    tdtFilt = filters.gaussian(tdt, sigma=sig2)\n",
    "    #plt.imshow(tdtFilt)\n",
    "    #plt.show()\n",
    "    #plt.imshow(dapiFilt)\n",
    "    #plt.show()\n",
    "    \n",
    "    #Apply mask\n",
    "    dapiMask = dapiFilt > cut1\n",
    "    tdtMask = tdtFilt > cut2\n",
    "    \n",
    "    #Plot if desired\n",
    "    if plot:\n",
    "        plt.imshow(tdtMask)\n",
    "        plt.show()\n",
    "        plt.imshow(dapiMask)\n",
    "        plt.show()\n",
    "    if plot:\n",
    "        flat_tdt = tdtFilt.flatten()\n",
    "        dist = plt.hist(flat_tdt[np.where(flat_tdt > 0)], density = True, color = \"magenta\", alpha = 0.5, bins = 1000, cumulative = False)\n",
    "        plt.show()\n",
    "    \n",
    "    #Compute the nestin profile\n",
    "    #Tdt+ regions are donor, Tdt- are host\n",
    "    nes_donor = (nes+1)*dapiMask*tdtMask\n",
    "    nes_host = (nes+1)*dapiMask*~tdtMask\n",
    "    \n",
    "    #Compute profile for gene of interest\n",
    "    goi_donor = (goi+1)*dapiMask*tdtMask\n",
    "    goi_host = (goi+1)*dapiMask*~tdtMask\n",
    "    \n",
    "    #Plot in the specified region\n",
    "    #If both are zero, then want for whole image\n",
    "    if r1 == 0 and r2 == 0:\n",
    "        nes_host_flat = nes_host.flatten()\n",
    "        nes_donor_flat = nes_donor.flatten()\n",
    "        goi_host_flat = goi_host.flatten()\n",
    "        goi_donor_flat = goi_donor.flatten()\n",
    "    \n",
    "    #If just r1 is zero, then it is just a stripe\n",
    "    elif r1 == 0:\n",
    "        nes_host_flat = nes_host[:, r2[0]:r2[1]].flatten()\n",
    "        nes_donor_flat = nes_donor[:, r2[0]:r2[1]].flatten()\n",
    "        goi_host_flat = goi_host[:, r2[0]:r2[1]].flatten()\n",
    "        goi_donor_flat = goi_donor[:, r2[0]:r2[1]].flatten()\n",
    "    \n",
    "    #Also stripe for r2\n",
    "    elif r2 == 0:\n",
    "        nes_host_flat = nes_host[r1[0]:r1[1], :].flatten()\n",
    "        nes_donor_flat = nes_donor[r1[0]:r1[1], :].flatten()\n",
    "        goi_host_flat = goi_host[r1[0]:r1[1], :].flatten()\n",
    "        goi_donor_flat = goi_donor[r1[0]:r1[1], :].flatten()\n",
    "        \n",
    "    #Otherwise, we want everything in the well-defined rectangle\n",
    "    else:\n",
    "        nes_host_flat = nes_host[r1[0]:r1[1], r2[0]:r2[1]].flatten()\n",
    "        nes_donor_flat = nes_donor[r1[0]:r1[1], r2[0]:r2[1]].flatten()\n",
    "        goi_host_flat = goi_host[r1[0]:r1[1], r2[0]:r2[1]].flatten()\n",
    "        goi_donor_flat = goi_donor[r1[0]:r1[1], r2[0]:r2[1]].flatten()\n",
    "    \n",
    "    #Get only non-zero pixels\n",
    "    nes_host_flat = nes_host_flat[np.where(nes_host_flat > 0)]\n",
    "    nes_donor_flat = nes_donor_flat[np.where(nes_donor_flat > 0)]\n",
    "    \n",
    "    goi_host_flat = goi_host_flat[np.where(goi_host_flat > 0)]\n",
    "    goi_donor_flat = goi_donor_flat[np.where(goi_donor_flat > 0)]\n",
    "    \n",
    "    #Flip the host and donor colors to get them to match the right species\n",
    "    if \"M-R\" in file:\n",
    "        host_color = \"green\"\n",
    "        donor_color = \"magenta\"\n",
    "    elif \"R-M\" in file:\n",
    "        host_color = \"magenta\"\n",
    "        donor_color = \"green\"\n",
    "    if plot:\n",
    "        #Plot histogram of expression for Nestin\n",
    "        plt.hist(np.log10(nes_donor_flat), density = True, color = donor_color, alpha = 0.5, bins = 100)\n",
    "        plt.hist(np.log10(nes_host_flat), density = True, color = host_color, alpha = 0.5, bins = 100)\n",
    "        plt.title(\"Nestin expression\")\n",
    "        plt.show()\n",
    "        \n",
    "        #Do the same for the gene of interest\n",
    "        plt.hist(np.log10(goi_donor_flat), density = True, color = donor_color, alpha = 0.5, bins = 100)\n",
    "        plt.hist(np.log10(goi_host_flat), density = True, color = host_color, alpha = 0.5, bins = 100)\n",
    "        plt.title(file.split(\"/\")[2].replace(\".czi\", \"\") + \" expression\")\n",
    "        plt.show()\n",
    "        \n",
    "        #Print the median value\n",
    "        print(np.median(goi_host_flat))\n",
    "    \n",
    "    #Compute background subtraction\n",
    "    goi_back = np.mean(goi[back1:, 0:back2].flatten() + 1)\n",
    "    \n",
    "    #Return background subtracted arrays for the gene of interest\n",
    "    return goi_host_flat-goi_back, goi_donor_flat-goi_back\n",
    "\n",
    "#Plot across all conditions\n",
    "#Takes as input the arrays returned by plot defined above\n",
    "def plot_across(mr_host, mr_donor, rm_host, rm_donor, gene):\n",
    "    \n",
    "    #Get only non-zero ones and log transform\n",
    "    #Then plot histogram of pixel values\n",
    "    plt.hist(np.log10(mr_host[np.where(mr_host > 1)]), density = True, color = \"green\", alpha = 0.5, bins = 100)\n",
    "    plt.hist(np.log10(rm_host[np.where(rm_host > 1)]), density = True, color = \"magenta\", alpha = 0.5, bins = 100)\n",
    "    plt.title(gene + \" expression hosts\")\n",
    "    plt.show()\n",
    "    plt.hist(np.log10(mr_donor[np.where(mr_donor > 1)]), density = True, color = \"magenta\", alpha = 0.5, bins = 100)\n",
    "    plt.hist(np.log10(rm_donor[np.where(rm_donor > 1)]), density = True, color = \"green\", alpha = 0.5, bins = 100)\n",
    "    plt.title(gene + \" expression donors\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d726f66f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Iterate through computing based on different masks\n",
    "#Designed to make sure that the estimates are not sensitive to the parameters chosen\n",
    "\n",
    "#Define masks for RM and MR images\n",
    "rm_dapi_masks = [0.08, 0.1, 0.12]\n",
    "rm_tdt_masks = [0.05, 0.06, 0.07]\n",
    "mr_dapi_masks = [0.1, 0.12, 0.14]\n",
    "mr_tdt_masks = [0.1, 0.11, 0.12]\n",
    "out = []\n",
    "\n",
    "#Iterate through all combinations\n",
    "for rm_dapi_mask in rm_dapi_masks:\n",
    "    for rm_tdt_mask in rm_tdt_masks:\n",
    "        for mr_dapi_mask in mr_dapi_masks:\n",
    "            for mr_tdt_mask in mr_tdt_masks:\n",
    "                \n",
    "                #Can change file name to the correct file\n",
    "                rm_host, rm_donor = plot(\"Batch3/R-M Ch8/c-Jun.czi\", 10, 5, rm_dapi_mask, rm_tdt_mask, [800, 1600], [1400, 3300], 2750, 1000, plot = False)\n",
    "                mr_host, mr_donor = plot(\"Batch3/M-R Ch1/c-Jun.czi\", 10, 5, mr_dapi_mask, mr_tdt_mask, [1000, 3000], [2750, 5000], 4000, 1000, plot = False)\n",
    "                \n",
    "                #Compute mean non-zero log2 transformed expression\n",
    "                mr_host_mean = np.mean(np.log2(mr_host[np.where(mr_host > 1)]))\n",
    "                mr_donor_mean = np.mean(np.log2(mr_donor[np.where(mr_donor > 1)]))\n",
    "                rm_host_mean = np.mean(np.log2(rm_host[np.where(rm_host > 1)]))\n",
    "                rm_donor_mean = np.mean(np.log2(rm_donor[np.where(rm_donor > 1)]))\n",
    "                \n",
    "                #Append\n",
    "                out.append([rm_dapi_mask, rm_tdt_mask, mr_dapi_mask, mr_tdt_mask, mr_donor_mean - mr_host_mean, rm_host_mean - rm_donor_mean, rm_host_mean - mr_host_mean, mr_donor_mean - rm_donor_mean])\n",
    "df = pd.DataFrame(out)\n",
    "#Make column names\n",
    "df.columns = [\"RM_DAPI_mask\", \"RM_TdT_mask\", \"MR_DAPI_mask\", \"MR_TdT_mask\", \"DM/HR\", \"HM/DR\", \"HM/HR\", \"DM/DR\"]\n",
    "\n",
    "#Write out, can change file name\n",
    "df.to_csv(\"Batch3_RM-Ch8_MR-Ch1_c-Jun.csv\", sep = \",\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56e6ef30",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Print the standard errors\n",
    "#Low standard error suggests that it is not sensitive to our choice of parameters\n",
    "v = pd.read_csv(\"Batch3_RM-Ch8_MR-Ch1_c-Jun.csv\", sep = \",\")\n",
    "print(np.sqrt(np.var(v[\"HM/HR\"]))/9)\n",
    "print(np.sqrt(np.var(v[\"DM/HR\"]))/9)\n",
    "print(np.sqrt(np.var(v[\"DM/DR\"]))/9)\n",
    "\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c2c9dae",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Actual parameters used for Jun RM\n",
    "rm_host, rm_donor = plot(\"Batch3/R-M Ch8/c-Jun.czi\", 10, 5, 0.1, 0.07, [800, 1600], [1400, 3300], 2750, 1000, plot = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "045b6ee1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Actual parameters used for Jun MR\n",
    "mr_host, mr_donor = plot(\"Batch3/M-R Ch1/c-Jun.czi\", 10, 5, 0.14, 0.12, [1000, 3000], [2750, 5000], 4000, 1000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51d15d2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot everything to check it is okay\n",
    "plot_across(mr_host, mr_donor, rm_host, rm_donor, \"c-Jun\")\n",
    "\n",
    "#Compute mean of log2 transformed non-zero data\n",
    "mr_host_mean = np.mean(np.log2(mr_host[np.where(mr_host > 1)]))\n",
    "mr_donor_mean = np.mean(np.log2(mr_donor[np.where(mr_donor > 1)]))\n",
    "rm_host_mean = np.mean(np.log2(rm_host[np.where(rm_host > 1)]))\n",
    "rm_donor_mean = np.mean(np.log2(rm_donor[np.where(rm_donor > 1)]))\n",
    "\n",
    "#Print the log2 fold-changes needed to compute what we want\n",
    "print(\"DM/HR\", mr_donor_mean - mr_host_mean)\n",
    "print(\"HM/DR\", rm_host_mean - rm_donor_mean)\n",
    "print(\"HM/HR\", rm_host_mean - mr_host_mean)\n",
    "print(\"DM/DR\", mr_donor_mean - rm_donor_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82d8f4a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Iterate through computing based on different masks,this time for batch1 of Jun measurements\n",
    "rm_dapi_masks = [0.09, 0.1, 0.11]\n",
    "rm_tdt_masks = [0.035, 0.04, 0.045]\n",
    "mr_dapi_masks = [0.16, 0.18, 0.20]\n",
    "mr_tdt_masks = [0.11, 0.13, 0.15]\n",
    "out = []\n",
    "for rm_dapi_mask in rm_dapi_masks:\n",
    "    for rm_tdt_mask in rm_tdt_masks:\n",
    "        for mr_dapi_mask in mr_dapi_masks:\n",
    "            for mr_tdt_mask in mr_tdt_masks:\n",
    "                rm_host, rm_donor = plot(\"Batch1/R-M Ch2/c-Jun.czi\", 10, 5, rm_dapi_mask, rm_tdt_mask, [800, 1600], [1400, 3300], 2750, 1000, plot = False)\n",
    "                mr_host, mr_donor = plot(\"Batch1/M-R Ch5/c-Jun.czi\", 10, 5, mr_dapi_mask, mr_tdt_mask, [1000, 3000], [2750, 5000], 4000, 1000, plot = False)\n",
    "                \n",
    "                mr_host_mean = np.mean(np.log2(mr_host[np.where(mr_host > 1)]))\n",
    "                mr_donor_mean = np.mean(np.log2(mr_donor[np.where(mr_donor > 1)]))\n",
    "                rm_host_mean = np.mean(np.log2(rm_host[np.where(rm_host > 1)]))\n",
    "                rm_donor_mean = np.mean(np.log2(rm_donor[np.where(rm_donor > 1)]))\n",
    "                \n",
    "                out.append([rm_dapi_mask, rm_tdt_mask, mr_dapi_mask, mr_tdt_mask, mr_donor_mean - mr_host_mean, rm_host_mean - rm_donor_mean, rm_host_mean - mr_host_mean, mr_donor_mean - rm_donor_mean])\n",
    "df = pd.DataFrame(out)\n",
    "df.columns = [\"RM_DAPI_mask\", \"RM_TdT_mask\", \"MR_DAPI_mask\", \"MR_TdT_mask\", \"DM/HR\", \"HM/DR\", \"HM/HR\", \"DM/DR\"]\n",
    "df.to_csv(\"Batch1_RM-Ch2_MR-Ch5_c-Jun.csv\", sep = \",\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00d6d16e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Print Standard deviations again\n",
    "v = pd.read_csv(\"Batch1_RM-Ch2_MR-Ch5_c-Jun.csv\", sep = \",\")\n",
    "print(np.sqrt(np.var(v[\"HM/HR\"]))/9)\n",
    "print(np.sqrt(np.var(v[\"DM/HR\"]))/9)\n",
    "print(np.sqrt(np.var(v[\"DM/DR\"]))/9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59d9f0c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Repeat above for batch1 values\n",
    "mr_host, mr_donor = plot(\"Batch1/M-R Ch5/c-Jun.czi\", 10, 5, 0.16, 0.11, [100, 1000], [4300, 5800], 3000, 1000)\n",
    "\n",
    "rm_host, rm_donor = plot(\"Batch1/R-M Ch2/c-Jun.czi\", 10, 5, 0.11, 0.04, [1300, 2250], [2500, 3500], 2000, 1000)\n",
    "\n",
    "mr_host_mean = np.mean(np.log2(mr_host[np.where(mr_host > 1)]))\n",
    "mr_donor_mean = np.mean(np.log2(mr_donor[np.where(mr_donor > 1)]))\n",
    "rm_host_mean = np.mean(np.log2(rm_host[np.where(rm_host > 1)]))\n",
    "rm_donor_mean = np.mean(np.log2(rm_donor[np.where(rm_donor > 1)]))\n",
    "\n",
    "print(mr_donor_mean - mr_host_mean)\n",
    "print(rm_host_mean - rm_donor_mean)\n",
    "print(rm_host_mean - mr_host_mean)\n",
    "print(mr_donor_mean - rm_donor_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef683a31",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now do computations for Hspa5 as done for Jun\n",
    "mr_host, mr_donor = plot(\"Batch1/M-R Ch5/Hspa5.czi\", 10, 5, 0.2, 0.15, [150, 1100], [4400, 5700], 3000, 1000)\n",
    "rm_host, rm_donor = plot(\"Batch1/R-M Ch2/Hspa5.czi\", 10, 5, 0.15, 0.032, [1300, 2250], [1250, 3500], 2000, 1000)\n",
    "\n",
    "plot_across(mr_host, mr_donor, rm_host, rm_donor, \"Hspa5\")\n",
    "mr_host_mean = np.mean(np.log2(mr_host[np.where(mr_host > 1)]))\n",
    "mr_donor_mean = np.mean(np.log2(mr_donor[np.where(mr_donor > 1)]))\n",
    "rm_host_mean = np.mean(np.log2(rm_host[np.where(rm_host > 1)]))\n",
    "rm_donor_mean = np.mean(np.log2(rm_donor[np.where(rm_donor > 1)]))\n",
    "\n",
    "print(\"DM/HR\", mr_donor_mean - mr_host_mean)\n",
    "print(\"HM/DR\", rm_host_mean - rm_donor_mean)\n",
    "print(\"HM/HR\", rm_host_mean - mr_host_mean)\n",
    "print(\"DM/DR\", mr_donor_mean - rm_donor_mean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e653ad51",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Can also similarly compute for Hspa5 Batch 3\n",
    "mr_host, mr_donor = plot(\"Batch3/M-R Ch1/Hspa5.czi\", 10, 5, 0.2, 0.2, [1000, 3000], [2750, 5000], 4000, 1000)\n",
    "rm_host, rm_donor = plot(\"Batch3/R-M Ch8/Hspa5.czi\", 10, 5, 0.1, 0.045, [750, 1500], [1250, 3750], 2750, 1000)\n",
    "plot_across(mr_host, mr_donor, rm_host, rm_donor, \"Hspa5\")\n",
    "mr_host_mean = np.mean(np.log2(mr_host[np.where(mr_host > 1)]))\n",
    "mr_donor_mean = np.mean(np.log2(mr_donor[np.where(mr_donor > 1)]))\n",
    "rm_host_mean = np.mean(np.log2(rm_host[np.where(rm_host > 1)]))\n",
    "rm_donor_mean = np.mean(np.log2(rm_donor[np.where(rm_donor > 1)]))\n",
    "\n",
    "print(\"DM/HR\", mr_donor_mean - mr_host_mean)\n",
    "print(\"HM/DR\", rm_host_mean - rm_donor_mean)\n",
    "print(\"HM/HR\", rm_host_mean - mr_host_mean)\n",
    "print(\"DM/DR\", mr_donor_mean - rm_donor_mean)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}