scripts/Figure_5_S6_Tre1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Figures 5, S6\n",
    "This script runs our analysis of Tre1 intoxication of *E. coli* by *S. proteamaculans*. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "!nvcc --version\n",
    "!nvidia-smi\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import numpy as np\n",
    "import time, os, sys\n",
    "from urllib.parse import urlparse\n",
    "from urllib.parse import urlparse\n",
    "from cellpose_omni import models, core, io\n",
    "from skimage import measure\n",
    "from sklearn.cluster import k_means\n",
    "\n",
    "use_GPU = core.use_gpu()\n",
    "print('>>> GPU activated? %d'%use_GPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time, os, sys\n",
    "import skimage.io\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('dark_background')\n",
    "import matplotlib as mpl\n",
    "%matplotlib inline\n",
    "mpl.rcParams['figure.dpi'] = 300\n",
    "from cellpose_omni import utils, io, transforms\n",
    "# from cupyx import scipy\n",
    "\n",
    "# np.abs(np.array([1+1j,1]))\n",
    "# np.asarray(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "home = str(Path.home())\n",
    "basedir = '/home/kcutler/DataDrive/ecVSsp/'\n",
    "phase = Path(basedir).rglob('*c1.tif')\n",
    "pfiles = [str(p) for p in phase]\n",
    "gfiles = [p[:-5]+'2.tif' for p in pfiles]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = pfiles\n",
    "imgs = [skimage.io.imread(f) for f in pfiles]\n",
    "gfps = [skimage.io.imread(f) for f in gfiles]\n",
    "nimg = len(imgs)\n",
    "print(nimg)\n",
    "from cellpose_omni import transforms\n",
    "from skimage.filters import gaussian\n",
    "# plt.figure(figsize=(10,10))\n",
    "\n",
    "def localnormalize(im,sigma):\n",
    "    return im/gaussian(im,sigma)\n",
    "\n",
    "\n",
    "for k in range(len(imgs)):\n",
    "    img = transforms.move_min_dim(imgs[k])\n",
    "    gfp = transforms.move_min_dim(gfps[k])\n",
    "    if len(img.shape)>2:\n",
    "        imgs[k] = img[:,:,0]\n",
    "        gfps[k] = gfp[:,:,0]\n",
    "    imgs[k] = transforms.normalize99(imgs[k],omni=True)\n",
    "    # gfps[k] = localnormalize(gfps[k],sigma=100)\n",
    "    gfps[k] = transforms.normalize99(gfps[k],omni=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check that everything matches\n",
    "i = 5\n",
    "g = gfps[i].copy()\n",
    "g = g-gaussian(g,15)\n",
    "\n",
    "fig = plt.figure(figsize=(20,20))\n",
    "plt.imshow(np.hstack((imgs[i],gfps[i],g)))\n",
    "# p = plt.hist(g[g>.01],bins=100)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CP corresponds to 0, omni to 1\n",
    "# modeldir = ['/home/kcutler/DataDrive/omnipose_train/registered/models/cellpose_residual_on_style_on_concatenation_off_registered_2021_07_22_11_48_00.465748',\n",
    "#             '/home/kcutler/DataDrive/omnipose_train/registered/models/cellpose_residual_on_style_on_concatenation_off_registered_2021_09_30_21_41_17.488619_epoch_3999']\n",
    "model_type = ['bact_phase','bact_phase_omni']\n",
    "model = [models.CellposeModel(gpu=use_GPU, model_type=model_type[k],nchan=2) for k in range(2)]\n",
    "suffix = ['_CP','_OP']\n",
    "names = ['Cellpose','Omnipose']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "imglist = imgs[:] #easier selection\n",
    "n = len(names)\n",
    "N = len(imgs)\n",
    "J = range(len(model))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "masks,flows,styles,d = [[]]*n,[[]]*n,[[]]*n,[[]]*n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "clean = False\n",
    "if clean:\n",
    "    chans = [0,0]\n",
    "    for j in J:\n",
    "        masks[j], flows[j], styles[j] = model[j].eval(imglist,channels=chans,\n",
    "                                                      rescale=None,mask_threshold=-1,transparency=True,\n",
    "                                                      flow_threshold=0,omni=j,resample=False,tile=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "if clean:\n",
    "    for j in J:\n",
    "        savedir = '/home/kcutler/DataDrive/omnipose_paper/Figure 5/'\n",
    "        if not os.path.isdir(savedir):\n",
    "            os.mkdir(savedir)\n",
    "        io.save_masks(imglist, masks[j], flows[j], pfiles, suffix = suffix[j],save_flows=1,\n",
    "                      save_outlines=1,savedir=savedir,in_folders=True,save_txt=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "## Load masks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Already run segmentation\n",
    "def getname(path,suffix=''):\n",
    "    return os.path.splitext(Path(path).name)[0].replace(suffix,'')\n",
    "pnames = [getname(p) for p in pfiles]\n",
    "\n",
    "#add a line to read back in masks here\n",
    "maskdir = '/home/kcutler/DataDrive/omnipose_paper/Figure 5/masks/'\n",
    "om_files = [str(s) for s in Path(maskdir).rglob('*masks_OP.png')]\n",
    "cp_files = [str(s) for s in Path(maskdir).rglob('*masks_CP.png')]\n",
    "om_names = [getname(s,suffix='_cp_masks_OP') for s in om_files]\n",
    "cp_names = [getname(s,suffix='_cp_masks_CP') for s in cp_files]\n",
    "om_ind = np.array([ [idx for idx,s in enumerate(om_names) if s==p] for p in pnames]).flatten()\n",
    "cp_ind = np.array([ [idx for idx,s in enumerate(cp_names) if s==p] for p in pnames]).flatten()\n",
    "om_files = [om_files[i] for i in om_ind]\n",
    "cp_files = [cp_files[i] for i in cp_ind]\n",
    "\n",
    "om_masks = [skimage.io.imread(f)for f in om_files]\n",
    "cp_masks = [skimage.io.imread(f)for f in cp_files]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "#stardist masks\n",
    "maskdir_sd = '/home/kcutler/DataDrive/omnipose_paper/Figure 5/StarDist/'\n",
    "sd_files = [str(s) for s in Path(maskdir_sd).rglob('*.tif')]\n",
    "sd_names = [getname(s,suffix='_stardist_masks') for s in sd_files]\n",
    "sd_ind = np.array([ [idx for idx,s in enumerate(sd_names) if s==p] for p in pnames]).flatten()\n",
    "sd_files = [sd_files[i] for i in sd_ind]\n",
    "sd_masks = [skimage.io.\n",
    "            imread(f) for f in sd_files]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {},
   "outputs": [],
   "source": [
    "#MiSiC masks\n",
    "maskdir_ms = '/home/kcutler/DataDrive/omnipose_paper/Figure 5/MiSiC/'\n",
    "ms_files = [str(s) for s in Path(maskdir_ms).rglob('*.tif')]\n",
    "ms_names = [getname(s,suffix='_masks') for s in ms_files]\n",
    "ms_ind = np.array([ [idx for idx,s in enumerate(ms_names) if s==p] for p in pnames]).flatten()\n",
    "ms_files = [ms_files[i] for i in ms_ind]\n",
    "ms_masks = [skimage.io.imread(f) for f in ms_files]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "## Define experimental categories \n",
    "Image set contains both initial time point and 20hr of wild-type and mutant (control) co-cultures. Here we separate the images into four groups, and only analyse two of them (20hr time point). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['StarDist','Cellpose','Omnipose','MiSiC']\n",
    "abbrev =  ['SD','CP','OP','MS']\n",
    "masks = [sd_masks,cp_masks,om_masks,ms_masks]\n",
    "J = range(len(masks))\n",
    "# plt.imshow(flows[0][0][0])\n",
    "def getname(path,suffix=''):\n",
    "        return os.path.splitext(Path(path).name)[0].replace(suffix,'')\n",
    "img_names = [getname(f) for f in pfiles]\n",
    "\n",
    "cat1 = ['wt','mut']\n",
    "cat2 = '20hr'\n",
    "indices = [[[] for j in range(2)] for k in range(2)]\n",
    "K = range(len(cat1))\n",
    "for k in K:\n",
    "    indices[k][0] = [i for i, s in enumerate(pfiles) if (cat1[k] in s) and (cat2 not in s)]\n",
    "    indices[k][1] = [i for i, s in enumerate(pfiles) if (cat1[k] in s) and (cat2 in s)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "pnames"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "## Collect data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "intensity = [[[] for k in K] for j in J]\n",
    "intensity_eroded = [[[] for k in K] for j in J]\n",
    "area = [[[] for k in K] for j in J]\n",
    "regions = [[[] for k in K] for j in J]\n",
    "species = [[[] for k in K] for j in J]\n",
    "\n",
    "from scipy.ndimage import binary_erosion\n",
    "from skimage.morphology import thin\n",
    "for j in J:\n",
    "    print(names[j])\n",
    "    for k in K:\n",
    "        inds = indices[k][1] #[1]->just process 20hr images\n",
    "        for i in inds: # loop over each group of images \n",
    "            regs = measure.regionprops(masks[j][i],intensity_image=gfps[i])\n",
    "            for r in regs[1:]: #no background\n",
    "                # mask = binary_erosion(r.image,iterations=2)\n",
    "                mask = thin(r.image,max_num_iter=6)\n",
    "                if mask is not None and np.sum(mask)>0: \n",
    "                    r.J = j\n",
    "                    r.K = k\n",
    "                    r.img_index = i\n",
    "                    r.img_name = img_names[i]\n",
    "                    intensity[j][k].append(r.mean_intensity)\n",
    "                    intensity_eroded[j][k].append(np.mean(r.intensity_image[mask])) #avoid edge artifacts from bleedover/GFP registration \n",
    "                    area[j][k].append(r.area)\n",
    "                    regions[j][k].append(r)\n",
    "                    if len(area[j][k])!=len(intensity[j][k]):\n",
    "                        print(i,j,k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(area[j][k])==len(intensity[j][k]) # sanity check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "text = ['Tre 1 experiment','image count: {}'.format(len(imgs)),'']\n",
    "\n",
    "for j in J:\n",
    "    text += ['{} masks\\n'.format(names[j])]\n",
    "    cellcounts = [len(area[j][k]) for k in K]\n",
    "    for k in K:\n",
    "        text+=['\\t{} cell count: {}'.format(cat1[k],cellcounts[k])]\n",
    "    text+=['\\ttotal cell count: {}\\n'.format(np.sum(cellcounts))]\n",
    "\n",
    "basedir = '/home/kcutler/DataDrive/omnipose_all/'        \n",
    "with open(os.path.join(basedir,'Tre1_stats.txt'), \"w\") as text_file:\n",
    "    [print(t,file=f) for t in text for f in [None,text_file]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "J = range(len(masks))\n",
    "I= [[] for j in J]\n",
    "A = [[] for j in J]\n",
    "R = [[] for j in J]\n",
    "pixelArea = 0.065**2\n",
    "for j in J:\n",
    "    # I[j] = [np.array([item for sublist in intensity_eroded[j][k] for item in sublist]) for k in K]\n",
    "    # A[j] = [np.array([item for sublist in area[j][k] for item in sublist])*pixelArea for k in K]\n",
    "    # R[j] = [np.array([item for sublist in regions[j][k] for item in sublist]) for k in K]\n",
    "    I[j] = [np.nan_to_num(intensity_eroded[j][k]) for k in K]\n",
    "    A[j] = [np.array(area[j][k]) for k in K]\n",
    "    R[j] = [np.array(regions[j][k]) for k in K]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# len(I[0][0]),len(I[1][0]),len(I[2][0])\n",
    "# np.argwhere(A[0][1]==A[1][1])\n",
    "bins=50\n",
    "# p = plt.hist(I[1],bins=bins) \n",
    "j = 2\n",
    "print(names[j])\n",
    "fig,ax = plt.subplots()\n",
    "p2 = plt.hist(I[j],bins=bins,density=True)\n",
    "ax.set_xscale('log')\n",
    "# len(intensity[1][1]),len(intensity[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Is = np.concatenate([I[j][k]*A[j][k] for j,k in zip(J,K)])\n",
    "# As = np.concatenate([A[j][k] for j,k in zip(J,K)])\n",
    "# mI = np.median(Is)/np.median(As)\n",
    "J = range(len(masks))\n",
    "Is = [I[j][k] for j,k in zip(J,K)]\n",
    "# Is = (Is-Is.min())/(Is.max()-Is.min())\n",
    "mI=[np.percentile(Is[k],50) for k in K]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define species per mask\n",
    "species = [[[] for k in K] for j in J]\n",
    "labels = [[[] for k in K] for j in J]\n",
    "        \n",
    "for j in J:    \n",
    "    for k in K:\n",
    "        i = I[j][k].copy()/mI[k]\n",
    "        X = i.reshape(-1, 1)\n",
    "        km = k_means(X,4,random_state=42) # k-means binning; km[1] is an index list of bins\n",
    "        idx = np.argsort(km[0].flatten())\n",
    "        # inds = np.logical_or(km[1]==idx[0], km[1]==idx[1])\n",
    "        inds = km[1]==idx[0]\n",
    "        species[j][k] = inds\n",
    "        labels[j][k] = [r.label for r in R[j][k]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {},
   "source": [
    "## Make Figure 5,S6 Segmentation plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cellpose_omni import plot, metrics\n",
    "import omnipose\n",
    "from omnipose.utils import rescale\n",
    "import ncolor\n",
    "from scipy.optimize import linear_sum_assignment\n",
    "import fastremap\n",
    "from scipy.ndimage import binary_dilation\n",
    "cmap = mpl.cm.get_cmap('plasma')\n",
    "cmap2 = mpl.cm.get_cmap('gray')\n",
    "outl_col = cmap(0.85)[:3]\n",
    "pad = 10\n",
    "bg = 0.5\n",
    "savedir = '/home/kcutler/DataDrive/omnipose_all/tre1/bigcells'\n",
    "io.check_dir(savedir)\n",
    "\n",
    "J = range(len(masks))\n",
    "# for k in K:\n",
    "#to plot these points on the graph, I either need to track a bunch of indices, or I make a suplemental list of areas \n",
    "#vs fluroescence (normalized also by mI). because I compressed R, I, and A, I don't think it will be easy to find\n",
    "#indices. So for each Omnipose long cell, we make a list of corresponding area/intesity points to plot. \n",
    "\n",
    "\n",
    "for k in [0]:\n",
    "    # big_ind = np.argwhere(A[j][k]>15).flatten()\n",
    "    \n",
    "    big_ind = [16550,29239,38647,51014,164489,203681,76074,100162,101979] #[5885,5901] #[5885,5901] 100162 \n",
    "    # big_ind = [big_ind[4]]#4,-2\n",
    "    big_ind = [big_ind[i] for i in [0,4,6]]#4,-2\n",
    "    dotcolors = ['#FF0800','#FF903E','#BBBDBF']\n",
    "    \n",
    "\n",
    "    # labels = [r.label for r in R[j][k][species[j][k]]]\n",
    "    labels0 = labels[j][k]\n",
    "    L = len(big_ind)\n",
    "    cell_error_list = [[[] for l in range(L)] for j in J]\n",
    "    for idx,i in enumerate(big_ind):\n",
    "    # for i in [big_ind[0]]:\n",
    "    # for i in [16557]:\n",
    "    # for i in range(16540,16557):\n",
    "        print(i)\n",
    "        j = 2 #omni \n",
    "        reg = R[j][k][i]\n",
    "        l = reg.label\n",
    "        print('Cell info',l,reg.area*pixelArea)\n",
    "        \n",
    "        # extract data\n",
    "        mask = masks[j][reg.img_index]\n",
    "        p = imgs[reg.img_index]\n",
    "        im = p**(np.log(bg)/np.log(np.mean(p[binary_erosion(mask==0)])))\n",
    "        gf = gfps[reg.img_index]\n",
    "        ly,lx = im.shape\n",
    "        \n",
    "        #put in coordinates of the cell of interest (Omnipose)\n",
    "        cell_error_list[2][idx] = [[A[j][k][i],I[j][k][i]]] #list with one\n",
    "        \n",
    "        #get single mask and all the neighbors too \n",
    "        cellmask = mask==l\n",
    "        touching = np.unique(mask[binary_dilation(cellmask)])\n",
    "        neighbor_labels = [ll for ll in labels0 if ll in touching]\n",
    "        neighbor_mask = np.zeros_like(cellmask) \n",
    "        for nl in neighbor_labels:\n",
    "            neighbor_mask[mask==nl] = nl\n",
    "            \n",
    "        y,x = np.nonzero(neighbor_mask)\n",
    "        y0 = max(np.min(y)-pad,0)\n",
    "        y1 = min(np.max(y)+pad,ly)\n",
    "        x0 = max(np.min(x)-pad,0)\n",
    "        x1 = min(np.max(x)+pad,lx)\n",
    "            \n",
    "        img0 = np.stack([im[y0:y1,x0:x1]]*3,axis=-1)\n",
    "        gfp0 = np.stack([gf[y0:y1,x0:x1]]*3,axis=-1)\n",
    "        om_m = mask[y0:y1,x0:x1] \n",
    "        # bin0 = np.logical_or(om_m==l,om_m==l-1) # capture both ong cells in this FoV\n",
    "        bin0 = om_m==l\n",
    "        a = np.count_nonzero(bin0)*pixelArea\n",
    "        name = 'area_'+str(round(a))+'_cell'+str(i)\n",
    "        outline = plot.outline_view(img0,om_m,color=outl_col,mode='thick')\n",
    "        skimage.io.imsave(os.path.join(savedir,name+'_'+abbrev[2]+'_outline'+'.png'),np.uint8(outline*255)) #omnipose\n",
    "        \n",
    "        mgt, remap = fastremap.renumber(om_m) #take omnipose as ground truth \n",
    "        # mgt = om_m.copy()\n",
    "        for ii in [1,0,-1]: #cellpose, stardist, misic\n",
    "            print(names[ii])\n",
    "            mpred, rmp = fastremap.renumber(masks[ii][reg.img_index][y0:y1,x0:x1])\n",
    "#             iou = metrics._intersection_over_union(mgt, mpred)\n",
    "#             th = 0\n",
    "#             n_min = min(iou.shape[0], iou.shape[1])\n",
    "#             costs = -(iou >= th).astype(float) - iou / (2*n_min)\n",
    "#             true_ind, pred_ind = linear_sum_assignment(costs)\n",
    "#             miou = iou[true_ind, pred_ind]\n",
    "            ind = remap[l]\n",
    "#             # ind = l\n",
    "#             # s_iou = iou[ind]\n",
    "#             # hits = np.argwhere(s_iou>.2)\n",
    "#             # argpart = np.argpartition(s_iou,3)[:3] # DON\"T DO IOU - need overlap code for errors\n",
    "    \n",
    "#             # print(ind,true_ind[ind],pred_ind[ind],miou[ind],hits)\n",
    "            \n",
    "            \n",
    "            fov_regs = measure.regionprops(mpred,intensity_image=gf[y0:y1,x0:x1])\n",
    "            ovp = metrics._label_overlap(mgt, mpred)[1:,1:] #throw out columns corresponding to zero  \n",
    "            areas = np.array([r.area for r in fov_regs])\n",
    "            # print(len(fov_regs),areas.shape,len(np.unique(mpred)),len(np.unique(mgt)),mgt.shape,mpred.shape,ovp.shape)\n",
    "            ovr = ovp / areas[np.newaxis,:] # Overlap Ratio\n",
    "            # p = plt.figure()\n",
    "            # plt.imshow(np.hstack((utils.rescale(ovp),utils.rescale(ovr>.5))))\n",
    "            # plt.imshow(utils.rescale(ovr))\n",
    "            # p.show()\n",
    "            # ind2 = np.argmax(np.sum(ovr>.75,axis=1))\n",
    "            s_ovr = ovr[ind-1]\n",
    "            # print(ind,ind2,s_ovr[s_ovr>0])#np.argmax(np.sum(ovr,axis=1))\n",
    "            hits = np.argwhere(s_ovr>.5).flatten() #these columns should correspond to the remapped prediction labels \n",
    "            # print(hits)\n",
    "            hitareas = np.array([areas[h] for h in hits]).flatten()\n",
    "            print('Areas',hitareas)\n",
    "\n",
    "\n",
    "            errors = []\n",
    "            hitmask = np.zeros_like(om_m) # for debugging\n",
    "            # for fov in fov_regs:\n",
    "            #     lb = fov.label\n",
    "            #     if lb in hits:\n",
    "            #         # print('sffd',lb in np.unique(mpred))\n",
    "            #         hitmask[mpred==lb+1] = lb+1 #label indexes are offset from real labels \n",
    "            #         a = fov.area\n",
    "            #         i = fov.mean_intensity\n",
    "            #         print('Point',a,i)\n",
    "            #         print('isinthere',a in hitareas)\n",
    "            #         errors.append([a*pixelArea,i/mI[k]])\n",
    "            \n",
    "            hit_regs = [r for r in fov_regs if r.area in hitareas] # label selection does not work; area does, but I worry about this (could over-select)\n",
    "            for fov in hit_regs:\n",
    "                \n",
    "                # sr,sc = reg.slice\n",
    "                # y0 = max(sr.start-pad,0)\n",
    "                # y1 = min(sr.stop+pad,ly)\n",
    "                # x0 = max(sc.start-pad-specialpad,0)\n",
    "                # x1 = min(sc.stop+pad,lx)\n",
    "                coords = fov.coords\n",
    "                y,x = zip(*(p for p in coords))\n",
    "                hitmask[y,x] = fov.label #label indexes are offset from real labels \n",
    "                a_ = fov.area\n",
    "                i_ = fov.mean_intensity\n",
    "                # print('Point',a,i)\n",
    "                # print('isinthere',a in hitareas)\n",
    "                errors.append([a_,i_])\n",
    "            # for lb in hits\n",
    "       \n",
    "            # for h in hits:\n",
    "            #     hitmask[mpred==h+1] = h+1\n",
    "\n",
    "            cell_error_list[ii][idx] = errors\n",
    "            \n",
    "            # if np.any(hitmask):\n",
    "            #     p = plt.figure(figsize=(10,10))\n",
    "            #     h = plot.outline_view(img0,hitmask)\n",
    "            #     ha = plot.outline_view(img0,mpred)\n",
    "            #     t = plot.outline_view(img0,mgt==ind)\n",
    "            #     plt.imshow(np.hstack((h,ha,t)))\n",
    "            #     plt.axis('off')\n",
    "            #     plt.show()\n",
    "            \n",
    "            \n",
    "            outline = plot.outline_view(img0,mpred,color=outl_col,mode='thick')\n",
    "            skimage.io.imsave(os.path.join(savedir,name+'_'+abbrev[ii]+'_outline'+'.png'),np.uint8(outline*255))\n",
    "            \n",
    "            \n",
    "\n",
    "        skimage.io.imsave(os.path.join(savedir,name+'_phase'+'.png'),np.uint8(img0*255))\n",
    "        skimage.io.imsave(os.path.join(savedir,name+'_gfp'+'.png'),np.uint8(gfp0*255))\n",
    "        \n",
    "        \n",
    "        ec_masks = [np.zeros_like(om_m)]*3\n",
    "        species_name = ['sp','ec']\n",
    "\n",
    "        if 1:\n",
    "            S = range(len(species_name))\n",
    "            species_masks = [[np.zeros_like(om_m) for j in J] for s in S] # want to separate Ec and Sp masks\n",
    "            for j in J:\n",
    "                print(names[j])\n",
    "                m = masks[j][reg.img_index][y0:y1,x0:x1]\n",
    "                reg0 = R[j][k][i]\n",
    "                # as of now, the plotting in the cell below needs to be run first for species partitioning...\n",
    "                #need to pull that out to run first\n",
    "                \n",
    "                labels1 = [[r.label for r in R[j][k][species[j][k]==s] if r.img_index == reg0.img_index and r.label in np.unique(m)] for s in S]\n",
    "                for s in S:\n",
    "                    for l in labels1[s]:\n",
    "                        species_masks[s][j][m==l] = 1\n",
    "                    \n",
    "                    nc = ncolor.label(m*species_masks[s][j])\n",
    "                    print(species_name[s])\n",
    "                    io.imsave(os.path.join(savedir,name+'_'+names[j]+'_mask_'+species_name[s]+'.png'),np.uint8(cmap2(1-rescale(nc))[:,:,:3]*255))\n",
    "                    plt.imshow(nc)\n",
    "                    plt.show()\n",
    "            # outline = plot.outline_view(img0,m,color=outl_col,mode='thick')\n",
    "            # io.imsave(os.path.join(savedir,name+'_'+names[j]+'_outline_'+species_name[s]+'.png'),np.uint8(outline*255))\n",
    "             \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "p = plt.figure()\n",
    "J = len(masks)\n",
    "for j in range(J):\n",
    "    errors = cell_error_list[j]\n",
    "    if errors:\n",
    "        for l in range(L):\n",
    "            pts = errors[l]\n",
    "            if pts: #if not empty\n",
    "                print(len(pts))\n",
    "                x, y = zip(*(p for p in pts))\n",
    "                plt.scatter(np.array(x)*pixelArea,y,facecolor=dotcolors[l],edgecolor=None)\n",
    "# cell_error_list[0]\n",
    "# p.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29",
   "metadata": {},
   "source": [
    "## Make Figure 5,S6 Scatter Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "labelsize = 9\n",
    "markersize = 5\n",
    "mpl.rcParams.update(mpl.rcParamsDefault)\n",
    "darkmode = 1\n",
    "if darkmode:\n",
    "    plt.style.use('dark_background')\n",
    "    axcol = 'k'\n",
    "    suffix = '_dark_mode'\n",
    "else:\n",
    "    mpl.rcParams.update(mpl.rcParamsDefault)\n",
    "    axcol = 'k'\n",
    "    # background_color = np.array([1,1,1,1])\n",
    "    suffix = ''\n",
    "\n",
    "background_color = [0,0,0,0]\n",
    "axcol = [0.5]*4\n",
    "    \n",
    "do_errors = False\n",
    "if not do_errors:\n",
    "    suffix+='_no_error_cell_colors'\n",
    "    \n",
    "J = len(masks)\n",
    "# sz = 3.25\n",
    "sz = 1\n",
    "xmax = 30\n",
    "# fig,axs = plt.subplots(J,2,figsize=(len(names)*sz,2*sz))\n",
    "# colors = [[['#8856a7','darkred'],['darkgray','dimgray']],\n",
    "#           [['limegreen','darkgreen'],['darkgray','dimgray']],\n",
    "#           [['orange','darkorange'],['darkgray','dimgray']],\n",
    "#          ]\n",
    "colors = [['#8BC53F','#58595B'],['#009345','#58595B']]\n",
    "# colors = ['#009345','#58595B']\n",
    "# cmap = mpl.cm.get_cmap('Greys')\n",
    "# cmap2 = mpl.cm.get_cmap('Greens')\n",
    "\n",
    "savedir = '/home/kcutler/DataDrive/omnipose_all/Tre1/plots/'\n",
    "\n",
    "io.check_dir(savedir)\n",
    "\n",
    "\n",
    "\n",
    "for j in range(J):\n",
    "    print(names[j])\n",
    "    errors = cell_error_list[j]\n",
    "\n",
    "    for k in K:\n",
    "# for j in [1]:\n",
    "#     for k in [1]:\n",
    "        # ax = axs[j][k]\n",
    "        fig,ax = plt.subplots(figsize=(3*sz,sz))\n",
    "        fig.patch.set_facecolor(background_color)\n",
    "        fig.patch.set_alpha(0)\n",
    "        ax.patch.set_facecolor(background_color)  # Set axes background to transparent\n",
    "        ax.set_rasterization_zorder(0)\n",
    "\n",
    "        a = A[j][k].copy()*pixelArea\n",
    "        i = I[j][k].copy()/mI[k]\n",
    "        # X = np.log(i).reshape(-1, 1)\n",
    "        X = i.reshape(-1, 1)\n",
    "        km = k_means(X,4,random_state=42,n_init=10) # k-means binning; km[1] is an index list of bins\n",
    "        idx = np.argsort(km[0].flatten())\n",
    "        # inds = np.logical_or(km[1]==idx[0], km[1]==idx[1])\n",
    "        inds = km[1]==idx[0]\n",
    "        species[j][k] = inds\n",
    "        ax.scatter(a[inds],i[inds],s=markersize,c=colors[k][1],rasterized=True, zorder=-1)\n",
    "        ax.scatter(a[~inds],i[~inds],s=markersize,c=colors[k][0],rasterized=True,zorder=-1)\n",
    "        \n",
    "        # if pts and k==0: #if not empty\n",
    "        #     x, y = zip(*(p for p in pts))\n",
    "        #     ax.scatter(x,y,c='r',s=markersize,rasterized=True)\n",
    "            \n",
    "        if errors and k==0 and do_errors:\n",
    "            for l in range(L):\n",
    "                pts = errors[l]\n",
    "                if pts: #if not empty\n",
    "                    x,y = zip(*(p for p in pts))\n",
    "                    a = np.array(x)*pixelArea\n",
    "                    i = np.array(y)/mI[k]\n",
    "                    ax.scatter(a,i,marker='o',facecolors=dotcolors[l], \n",
    "                               edgecolors='none',s=markersize*3,rasterized=False)\n",
    "#         scatter_heat(ax,a[inds],i[inds],s=markersize,cmap=cmap)\n",
    "#         scatter_heat(ax,a[~inds],i[~inds],s=markersize,cmap=cmap2)\n",
    "#         density_scatter(a[~inds],i[~inds],ax,s=markersize,cmap=cmap2)\n",
    "#         density_scatter(a[inds],i[inds],ax,s=markersize,cmap=cmap)\n",
    "\n",
    "        ax.tick_params(axis='both', which='major', labelsize=labelsize,\n",
    "                       length=3, direction=\"out\",colors=axcol,bottom=True,left=True)\n",
    "#         ax.set_xscale('log')\n",
    "        ax.set_yscale('log')\n",
    "        ax.set_ylim(.01,10)\n",
    "        ax.set_xlim(0,xmax) #110 covers all\n",
    "        ax.tick_params(which='both', colors=axcol)\n",
    "        ax.spines['right'].set_visible(False)\n",
    "        ax.spines['left'].set_color(axcol)\n",
    "        ax.spines['top'].set_visible(False)\n",
    "        ax.spines['bottom'].set_color(axcol)\n",
    "        \n",
    "        ax.patch.set_alpha(0.0)\n",
    "        # ax.set_title(names[j] + ' ' + cat1[k]+' '+str(a.size)+' %.2f'%(np.sum(inds)/a.size*100)) #+str(j)+str(k)\n",
    "        # ax.set_title(names[j] + ' ' + cat1[k] ) #+str(j)+str(k)\n",
    "        \n",
    "#         ax.set_axis_off()\n",
    "#         fig.add_axes(ax)\n",
    "        plt.show()\n",
    "        a = 0\n",
    "        tight_bbox_raw = ax.get_tightbbox(fig.canvas.get_renderer())\n",
    "        tight_bbox_raw._points+=[[-a,-a],[a,a]]\n",
    "        tight_bbox = mpl.transforms.TransformedBbox(tight_bbox_raw, mpl.transforms.Affine2D().scale(1./fig.dpi))\n",
    "        fig.savefig(os.path.join(savedir,names[j]+'_'+cat1[k]+suffix+'.pdf'),bbox_inches=tight_bbox, \n",
    "                    transparent=True, dpi=1000)\n",
    "#         ax.set_axis_on()\n",
    "# plt.subplots_adjust(hspace = .7,wspace = 0.25)\n",
    "# plt.show()\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31",
   "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}