dark_time_on_ring.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bbd216e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import os.path as ospath\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import shutil\n",
    "import glob\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "from picasso.picasso import io\n",
    "from picasso.picasso.postprocess import link, compute_dark_times\n",
    "from picasso.picasso.gui.render import estimate_kinetic_rate, fit_cum_exp\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4512a0dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_ring_data_df(dirname, filename):\n",
    "    file = ospath.join(dirname, filename)\n",
    "\n",
    "    try: \n",
    "        df = pd.read_pickle(file)\n",
    "    except FileNotFoundError:\n",
    "        print(\"The file {} does not exist.\".format(filename))\n",
    "        return None\n",
    "    else: \n",
    "        print(\"The file {} was loaded.\".format(filename))\n",
    "        return df\n",
    "    \n",
    "def load_data(path):\n",
    "\n",
    "    try:\n",
    "        locs, info = io.load_locs(path)\n",
    "    except io.NoMetadataFileError:\n",
    "        return None, None, None\n",
    "    \n",
    "    try:\n",
    "        pixelsize = info[1][\"Pixelsize\"]\n",
    "    except:  \n",
    "        print(\"No pixelsize found in yaml file. Default 130 nm used.\")\n",
    "        \n",
    "\n",
    "    if hasattr(locs, \"x_pick_rot\"):\n",
    "\n",
    "        # convert px to nm\n",
    "        locs.x *= pixelsize\n",
    "        locs.y *= pixelsize\n",
    "        locs.x_pick_rot *= pixelsize\n",
    "        locs.y_pick_rot *= pixelsize\n",
    "        \n",
    "        return locs, info, pixelsize\n",
    "    else:\n",
    "        print('x_pick_rot column is missing!')\n",
    "        return None, None, None\n",
    "\n",
    "        \n",
    "    #return locs, info, pixelsize\n",
    "\n",
    "\n",
    "def plot_cum_exp(pooled_locs, fit_result_len, fit_result_dark, locs_description, ax=None):\n",
    "    \n",
    "    color = blue\n",
    "\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    \n",
    "    data = pooled_locs.dark\n",
    "    data.sort()\n",
    "    y = np.arange(1, len(data) + 1)\n",
    "       \n",
    "    a = fit_result_dark.best_values[\"a\"]\n",
    "    t = fit_result_dark.best_values[\"t\"]\n",
    "    c = fit_result_dark.best_values[\"c\"]\n",
    "\n",
    "    ax.set_title(\n",
    "        \"Dark time (cumulative) - {} \\n\"\n",
    "        r\"$Fit: {:.2f}\\cdot(1-exp(x/{:.2f}))+{:.2f}$\".format(locs_description, a, t, c),loc=\"left\",fontsize=14)\n",
    "    data = pooled_locs.dark\n",
    "    data.sort()\n",
    "    y = np.arange(1, len(data) + 1)\n",
    "\n",
    "    ax.semilogx(data, y, c=color, label=\"Data\")\n",
    "    ax.semilogx(data, fit_result_dark.best_fit, c=red, label=\"Fit\")\n",
    "    ax.legend(loc=\"best\")\n",
    "    ax.set_xlabel(\"Duration (frames)\")\n",
    "    ax.set_ylabel(\"Frequency\")\n",
    "    #ax.set_xlim(1,data.max()+1)\n",
    "    \n",
    "    return ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1ba70976",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The file ring_data_filter.pkl was loaded.\n",
      "Index(['fov_id', 'cell_type', 'filename', 'group', 'residual', 'radius',\n",
      "       'angle_between_ring_and_coverslip', 'Sigma', 'n_events', 'mean_bright',\n",
      "       'mean_dark', 'fit_bright', 'fit_dark', 'filter_passed'],\n",
      "      dtype='object')\n"
     ]
    }
   ],
   "source": [
    "path = r'X:\\users\\kcramer\\sporPAINT\\ZapA\\ZapA_picked_rings'\n",
    "\n",
    "# consider locs within +- dR around the fitted ring radius\n",
    "dR = 100\n",
    "\n",
    "# linking of locs\n",
    "max_dist = 130 #nanometer\n",
    "max_dark_time = 15 #frames\n",
    "\n",
    "\n",
    "\n",
    "# load ring fit results\n",
    "df_ring_data = load_ring_data_df(path, \"ring_data_filter.pkl\")\n",
    "filtered = True\n",
    "if df_ring_data is None:\n",
    "    df_ring_data = load_ring_data_df(path, \"ring_data.pkl\")\n",
    "    filtered = False\n",
    "print(df_ring_data.keys())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "eb646cd9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1495786119f74231881d3b5a2076324a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Processing rings:   0%|          | 0/308 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "orientation = 'yz' # use yz, but also xy files contain info for filtering\n",
    "\n",
    "analysis_folder = os.path.join(path, \"analysis\")\n",
    "excluded_folder = os.path.join(path, 'analysis', 'excluded_filter')\n",
    "\n",
    "# image export settings\n",
    "img_format = \".png\"\n",
    "dpi = 100\n",
    "\n",
    "# define colors:\n",
    "blue = \"#4C72B0\"\n",
    "orange = \"#DD8452\"\n",
    "red = \"#C44E52\"\n",
    "gray = \"#90A8CE\"\n",
    "\n",
    "# prepare angle coordinates for plotting of rings check:\n",
    "# angles from -180 to 180\n",
    "bin_width = 360/12\n",
    "bins = np.arange(-180, 180+bin_width, bin_width)\n",
    "\n",
    "ring_data = []\n",
    "# add filter_passed column:\n",
    "df_ring_data['filename_ring'] =  'fov_' + df_ring_data['fov_id'].astype(str) +'_' + df_ring_data['cell_type'].astype(str) +'_pick_' + df_ring_data['group'].astype(str) +'_ring_' + '0' +'_rot_' + orientation + '_update.hdf5'\n",
    "for index, ring in tqdm(df_ring_data.iterrows(), desc=\"Processing rings\", total = df_ring_data.shape[0]):\n",
    "    plt.close('all')\n",
    "    locs, info, pixelsize = load_data(os.path.join(path, 'analysis', 'ring_locs', ring['filename_ring'] ))\n",
    "    \n",
    "    radius_fit = ring['radius']\n",
    "    radius_min = radius_fit - dR\n",
    "    radius_max = radius_fit + dR\n",
    "    \n",
    "    \n",
    "    N_locs = len(locs)\n",
    "\n",
    "    ################\n",
    "    # Get locs on ring\n",
    "    ################\n",
    "    locs_ring = locs[(locs['radius']>radius_min) & (locs['radius']<radius_max)]\n",
    "    df_locs_ring = pd.DataFrame.from_records(locs_ring)\n",
    "    locs_not_ring = locs[(locs['radius']<=radius_min) | (locs['radius']>=radius_max)]\n",
    "    df_locs_not_ring = pd.DataFrame.from_records(locs_not_ring)\n",
    "    \n",
    "\n",
    "\n",
    "    ################\n",
    "    # Link locs & calculate dark and bright times for all locs\n",
    "    ################\n",
    "    # link the isolated locs\n",
    "    r_locs_linked_all = link(locs, info, max_dist, max_dark_time)\n",
    "\n",
    "    # extract number of events\n",
    "    ring_n_events_all = len(r_locs_linked_all)\n",
    "    \n",
    "    # compute kinetics\n",
    "    kinetics_all = compute_dark_times(r_locs_linked_all)\n",
    "    #print('kinetics', kinetics)\n",
    "\n",
    "    # compute mean bright and dark times\n",
    "    ring_kinetics_all = {\"bright\":np.nanmean(kinetics_all.len),\n",
    "                         \"dark\":np.nanmean(kinetics_all.dark)}\n",
    "\n",
    "    fit_result_len_all = fit_cum_exp(kinetics_all.len)\n",
    "    fit_result_dark_all = fit_cum_exp(kinetics_all.dark)\n",
    "    \n",
    "    ################\n",
    "    # Link locs & calculate dark and bright times for locs on ring\n",
    "    ################\n",
    "    # link the isolated locs\n",
    "    r_locs_linked_ring = link(locs_ring, info, max_dist, max_dark_time)\n",
    "\n",
    "    # extract number of events\n",
    "    ring_n_events_ring = len(r_locs_linked_ring)\n",
    "    \n",
    "    # compute kinetics\n",
    "    kinetics_ring = compute_dark_times(r_locs_linked_ring)\n",
    "    #print('kinetics', kinetics)\n",
    "\n",
    "    # compute mean bright and dark times\n",
    "    ring_kinetics_ring = {\"bright\":np.nanmean(kinetics_ring.len),\n",
    "                         \"dark\":np.nanmean(kinetics_ring.dark)}\n",
    "\n",
    "    fit_result_len_ring = fit_cum_exp(kinetics_ring.len)\n",
    "    fit_result_dark_ring = fit_cum_exp(kinetics_ring.dark)\n",
    "    \n",
    "    \n",
    "    ################\n",
    "    # (5.2) append data into large array\n",
    "    ################\n",
    "\n",
    "    ring_kinetics_fit_all = {\"bright\":np.nanmean(fit_result_len_all.best_values[\"t\"]),\n",
    "                             \"dark\":np.nanmean(fit_result_dark_all.best_values[\"t\"])}\n",
    "    ring_kinetics_fit_ring = {\"bright\":np.nanmean(fit_result_len_ring.best_values[\"t\"]),\n",
    "                             \"dark\":np.nanmean(fit_result_dark_ring.best_values[\"t\"])}\n",
    "    \n",
    "    new_ring_data = [ring['fov_id'], #running file index\n",
    "                      ring['cell_type'], # veg or spore\n",
    "                      ring['filename'], #filename\n",
    "                      ring['group'], #pick number\n",
    "                      ring['residual'], #residual of ring fit\n",
    "                      ring['radius'], #fitted radius of spore (nm)\n",
    "                      ring['angle_between_ring_and_coverslip'], #anle between spore ring and coverslip (°)\n",
    "                      ring['Sigma'], #FWHM\n",
    "                      ring_n_events_all, #number of locs in spore ring\n",
    "                      ring_kinetics_all[\"bright\"], #mean bright time of spore ring (frames)\n",
    "                      ring_kinetics_all[\"dark\"], #mean dark time of spore ring (frames)\n",
    "                      ring_kinetics_fit_all[\"bright\"], #cumulative mean bright times of spore (frames)\n",
    "                      ring_kinetics_fit_all[\"dark\"], #cumulative mean dark times of spore (frames)\n",
    "                      ring_n_events_ring, #number of locs in spore ring\n",
    "                      ring_kinetics_ring[\"bright\"], #mean bright time of spore ring (frames)\n",
    "                      ring_kinetics_ring[\"dark\"], #mean dark time of spore ring (frames)\n",
    "                      ring_kinetics_fit_ring[\"bright\"], #cumulative mean bright times of spore (frames)\n",
    "                      ring_kinetics_fit_ring[\"dark\"] #cumulative mean dark times of spore (frames)\n",
    "                     ]\n",
    "    if filtered:\n",
    "        new_ring_data.append(ring['filter_passed'])\n",
    "        \n",
    "    ring_data.append(new_ring_data)\n",
    "    \n",
    "    \n",
    "    ################\n",
    "    # Plot\n",
    "    ################\n",
    "    \n",
    "    \n",
    "    # set up plot illustrating ring +- dR\n",
    "    fig = plt.figure(figsize=(20, 10), constrained_layout=True)\n",
    "    gs = fig.add_gridspec(1,3)\n",
    "    fig.suptitle((\"FOV {}, {}, Pick {} - $1/\\tau _d$ and $1/\\tau _b$ for locs on ring \\n\"\n",
    "                  \"File: {}\\n\").format(ring['fov_id'], ring['cell_type'] , ring['group'], ring['filename']), \n",
    "                 fontsize=16, \n",
    "                 ha=\"center\")\n",
    "\n",
    "    # set maximum radius for plot such that 99% of locs are shown\n",
    "    cutoff = 1\n",
    "\n",
    "    radius_sort = np.sort(locs['radius'])\n",
    "    radius_limit = radius_sort[int(cutoff*len(radius_sort))-1]\n",
    "\n",
    "    ## ring criterium\n",
    "    ax_1 = fig.add_subplot(gs[0, 0], projection='polar')\n",
    "    \n",
    "    ax_1.set_xticks(bins[1:]/180*np.pi)\n",
    "    ax_1.set_thetalim(-np.pi,np.pi)\n",
    "\n",
    "    #ring of fit radius\n",
    "    coords_angle_fit = np.arange(-np.pi,np.pi,0.01)\n",
    "    coords_radius_fit = np.full(len(coords_angle_fit), radius_fit)\n",
    "    coords_radius_min = np.full(len(coords_angle_fit), radius_min)\n",
    "    coords_radius_max = np.full(len(coords_angle_fit), radius_max)\n",
    "    \n",
    "    ax_1.fill_between(coords_angle_fit,radius_min, radius_max, color='red', alpha=0.2)\n",
    "\n",
    "    ax_1.plot(coords_angle_fit, coords_radius_fit, c='red', lw=2)\n",
    "    ax_1.plot(coords_angle_fit, coords_radius_min, c='red', lw=1)\n",
    "    ax_1.plot(coords_angle_fit, coords_radius_max, c='red', lw=1)\n",
    "    \n",
    "    # locs\n",
    "    #ax_1.scatter(locs['angle']/180*np.pi, locs['radius'], alpha=0.5)\n",
    "        # locs\n",
    "    ax_1.scatter(df_locs_not_ring['angle']/180*np.pi, df_locs_not_ring['radius'], alpha=0.1, c='#1f77b4')\n",
    "    ax_1.scatter(df_locs_ring['angle']/180*np.pi, df_locs_ring['radius'], alpha=0.5, c='#1f77b4')\n",
    "    \n",
    "    # set maximum radius for plot such that 99% of locs are shown\n",
    "    ax_1.set_rmax(radius_limit)\n",
    "\n",
    "    ax_1.set_title('Localizations within ring', fontsize = 16)\n",
    "    \n",
    "    \"\"\"\n",
    "    annotate_y = 0.85\n",
    "    plt.annotate(\"Fit radius: {:.1f} nm, defined ring width: $\\pm${:.0f} nm\\n\".expandtabs().format(radius_fit,dR),xy=(0.05,annotate_y-0.02), xycoords = 'subfigure fraction',fontsize=14)\n",
    "    plt.annotate(\"Percent of locs on ring: {:.1f}% \\t\\t Min. percent of locs on ring: {:.1f}%\\n\".expandtabs().format(perc_on_ring*100,perc_on_ring_min*100),xy=(0.05,annotate_y-0.04), xycoords = 'subfigure fraction',fontsize=14)\n",
    "    plt.annotate(\"Percent of locs inside of ring: {:.1f}%  \\t Max. percent of locs inside of ring: {:.1f}%\\n\".expandtabs().format(perc_inside_ring*100,perc_inside_ring_max*100),xy=(0.05,annotate_y-0.06), xycoords = 'subfigure fraction',fontsize=14)\n",
    "    plt.annotate(\"Percent of locs outside of ring: {:.1f}% \\t Max. percent of locs outside of ring: {:.1f}%\\n\".expandtabs().format(perc_outside_ring*100,perc_outside_ring_max*100),xy=(0.05,annotate_y-0.08), xycoords = 'subfigure fraction',fontsize=14)\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    ax2 = fig.add_subplot(gs[0, 1])\n",
    "    \n",
    "    plot_cum_exp(kinetics_all, fit_result_len_all, fit_result_dark_all,'all locs', ax=ax2)\n",
    "    \n",
    "    ax3 = fig.add_subplot(gs[0, 2])\n",
    "    \n",
    "    plot_cum_exp(kinetics_ring, fit_result_len_ring, fit_result_dark_ring, 'locs on ring', ax=ax3)\n",
    "    \n",
    "    \n",
    "    \n",
    "    img_fname = \"fov_{}_{}_pick_{}_kinetics_on_ring\".format(ring['fov_id'], ring['cell_type'], ring['group'])\n",
    "\n",
    "    plt.close(fig)\n",
    "    if filtered:\n",
    "        if ring['filter_passed'] == 'Yes':\n",
    "            img_name = os.path.join(analysis_folder, img_fname)\n",
    "            fig.savefig(img_name+img_format, dpi=dpi, format=\"png\")\n",
    "        elif ring['filter_passed'] == 'No':\n",
    "            img_name = os.path.join(excluded_folder, img_fname)\n",
    "            fig.savefig(img_name+img_format, dpi=dpi, format=\"png\")\n",
    "    else:\n",
    "        img_name = os.path.join(analysis_folder, img_fname)\n",
    "        fig.savefig(img_name+img_format, dpi=dpi, format=\"png\")\n",
    "        \n",
    "    \n",
    "\n",
    "df_ring_data = df_ring_data.drop(columns='filename_ring')\n",
    "\n",
    "\n",
    "columns = [\"fov_id\",\n",
    "            \"cell_type\",\n",
    "            \"filename\",\n",
    "            \"group\",\n",
    "            \"residual\",\n",
    "            \"radius\",\n",
    "            \"angle_between_ring_and_coverslip\",\n",
    "            \"Sigma\",\n",
    "            \"n_events_all\",\n",
    "            \"mean_bright_all\",\n",
    "            \"mean_dark_all\",\n",
    "            \"fit_bright_all\",\n",
    "            \"fit_dark_all\",\n",
    "            \"n_events_ring\",\n",
    "            \"mean_bright_ring\",\n",
    "            \"mean_dark_ring\",\n",
    "            \"fit_bright_ring\",\n",
    "            \"fit_dark_ring\",\n",
    "           ]\n",
    "\n",
    "if filtered:\n",
    "    columns.append('filter_passed')\n",
    "    \n",
    "# Save ring data\n",
    "df_ring_data = pd.DataFrame(ring_data, columns=columns)\n",
    "\n",
    "if filtered:\n",
    "    df_ring_data.to_csv(os.path.join(path,\"ring_data_filter_kinetics_on_ring.csv\"))\n",
    "    # Save dataframe with cell means for easy loading of data for postprocessing\n",
    "    df_ring_data.to_pickle(os.path.join(path,\"ring_data_filter_kinetics_on_ring.pkl\"))\n",
    "else:\n",
    "    df_ring_data.to_csv(os.path.join(path,\"ring_data_kinetics_on_ring.csv\"))\n",
    "    # Save dataframe with cell means for easy loading of data for postprocessing\n",
    "    df_ring_data.to_pickle(os.path.join(path,\"ring_data_kinetics_on_ring.pkl\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "440968bb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b86e00a3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}