In [1]:
# import the usual modules
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# nice plots
import seaborn as sns
sns.set(style="ticks") # turns the grid lines off 

# show images in viridis by default
plt.rcParams["image.cmap"] = "viridis"

# for getting file names
import glob

# skimage submodules we need
import skimage.io
import skimage.measure
import skimage.filters
import skimage.exposure

Part 1: using phase image to get segmentation mask

1.1 Background subtraction

In [2]:
# read in a sample image
im = skimage.io.imread("./data/lacI_titration/O2_R124_phase_pos_02.tif")

plt.imshow(im)
Out[2]:
<matplotlib.image.AxesImage at 0x1c17be9a20>
In [3]:
# normalize the image
im_norm = (im - np.min(im))/(np.max(im) - np.min(im))

# find the background image
gauss_radius = 50
im_bg = skimage.filters.gaussian(im_norm, gauss_radius)

# show this background image
plt.imshow(im_bg)
Out[3]:
<matplotlib.image.AxesImage at 0x1c189b49b0>
In [4]:
# perform background subtraction
im_bs = im_norm - im_bg

# show original the original and background subtracted side-by-side
fig, ax = plt.subplots(1,2)
ax[0].imshow(im_norm)
ax[1].imshow(im_bs)
Out[4]:
<matplotlib.image.AxesImage at 0x1c19301208>

1.2 Thresholding

Let's plot pixel intensities

In [5]:
hist, pixel_vals = skimage.exposure.histogram(im_bs)

plt.plot(pixel_vals, hist, '.')
plt.yscale("log")
plt.xlabel("pixel values")
plt.ylabel("counts")
Out[5]:
Text(0,0.5,'counts')
In [6]:
# choose a threshold based on plot above
threshold = -0.2

# threshold the image
im_thresh = im_bs < threshold

# show the thresholded image, bacteria are 1, backround is 0
plt.imshow(im_thresh)
Out[6]:
<matplotlib.image.AxesImage at 0x1c19d402e8>

1.3 label objects in the image and get their sizes

In [7]:
# label the objects in the binary image
im_label = skimage.measure.label(im_thresh)

# show the labeled image 
plt.imshow(im_label)
Out[7]:
<matplotlib.image.AxesImage at 0x1c19d9d780>
In [8]:
# find the number of object in the image
num_objs = np.max(im_label)
num_objs
Out[8]:
179

plot the size of the objects in the image

In [9]:
# length calibration
pixel_size = 0.16 # um

# array to store are values in 
areas = np.zeros(num_objs+1)

# obtain the areas of all labeled objects
for i in range(1,num_objs+1):
    # determine area for given object
    area = np.sum(im_label == i)
    areas[i] = area * pixel_size**2
In [10]:
# plot these areas
_ = plt.hist(areas, bins=20)
plt.yscale("log")
plt.xlabel("area (um)^2")
plt.ylabel("count")
Out[10]:
Text(0,0.5,'count')

1.4 Remove non-bacteria objects from the image

In [11]:
# set threshold on what we call a bacterium
lower_area = 1 # um^2
upper_area = 4 # um^2

# create a blank image with dimensions of im_label
im_mask = np.zeros_like(im_label)

# loop through objects of im_label
# add to im_mask *only* if the right size
for i in range(1,num_objs+1):
    # get area for given object
    pixel_area = np.sum(im_label == i)
    um_area = pixel_area * pixel_size**2
    
    # if area is within bounds, add this object to the mask
    if (um_area > lower_area) and (um_area < upper_area):
        im_mask = im_mask + (im_label == i)
In [12]:
# show the masked image
plt.imshow(im_mask)
Out[12]:
<matplotlib.image.AxesImage at 0x1c1a1c44a8>

1.5 "write" a masking function

In [13]:
import pboc_utils as pboc
In [14]:
test_mask = pboc.create_mask(im, area_high=4)

plt.imshow(test_mask)
Out[14]:
<matplotlib.image.AxesImage at 0x1c1a71b0b8>

Step 2. Use our segmentation mask to get fluorescence intensities from the corresponding YFP image

Let's show our phase, yfp, and mask image side by side

In [15]:
# read in the matching YFP image
im_yfp = skimage.io.imread("./data/lacI_titration/O2_R124_yfp_pos_02.tif")

# show images next to each other
fig, ax = plt.subplots(1,3, figsize=[10,3])
ax[0].imshow(im)

ax[1].imshow(im_yfp)

ax[2].imshow(test_mask)
Out[15]:
<matplotlib.image.AxesImage at 0x1c1b7d7828>

now let's determine the intensities of all the cells in the image

In [16]:
# list to store yfp intensities 
intensities = [] 

# number of cells in the mask
n_cells = np.max(test_mask)

# loop through cells and find the mean intensity for each
for i in range(1,n_cells):
    cell_mask = (test_mask == i)
    num_pixels = np.sum(cell_mask)
    cell_yfp = im_yfp * cell_mask
    cell_mean_intensity = np.sum(cell_yfp) / num_pixels
    intensities.append(cell_mean_intensity)
In [17]:
sns.swarmplot(intensities)
plt.xlabel("YFP expression (a.u.)")
Out[17]:
Text(0.5,0,'YFP expression (a.u.)')

example of how to write a function

In [18]:
def plus_one(x):
    ans = x + 1
    return ans

write a function to find the intensities of a given phase/yfp image pair

In [19]:
def find_intensities(im_phase,im_yfp):
    
    mask = pboc.create_mask(im_phase)
    
    # list to store yfp intensities 
    intensities = [] 

    # number of cells in the mask
    n_cells = np.max(mask)

    # loop through cells and find the mean intensity for each
    for i in range(1,n_cells):
        cell_mask = (mask == i)
        num_pixels = np.sum(cell_mask)
        cell_yfp = im_yfp * cell_mask
        cell_mean_intensity = np.sum(cell_yfp) / num_pixels
        intensities.append(cell_mean_intensity)
        
    return intensities

test out our function

In [20]:
find_intensities(im, im_yfp)
Out[20]:
[219.09174311926606,
 190.1818181818182,
 297.1980198019802,
 196.32,
 178.13483146067415,
 201.984496124031,
 194.183908045977,
 208.57943925233644,
 202.95145631067962,
 307.18674698795184,
 361.94871794871796,
 216.84057971014494,
 200.6904761904762,
 304.0263157894737,
 327.9313725490196,
 240.4121212121212]

Get the YFP intensities for all images of given strain (binding site and repressor copy number)

In [21]:
# strain information
operator = "O2"
repressor = "R124"

# filename structure
filename_phase = "./data/lacI_titration/" + operator + "_" + repressor + \
                "_phase_pos_*.tif"
filename_yfp = "./data/lacI_titration/" + operator + "_" + repressor + \
                "_yfp_pos_*.tif"

# grab all the images for this given strain
phase_names = np.sort(glob.glob(filename_phase))
yfp_names = np.sort(glob.glob(filename_yfp))

# number of positions
n_images = len(phase_names)

# list to store intensity values for *all* the images
intensities_all = []

# loop through the images
for i in range(n_images):
    
    # read in the image
    im_phase = skimage.io.imread(phase_names[i])
    im_yfp = skimage.io.imread(yfp_names[i])
    
    # get the intensities from from the YFP image
    intensities = find_intensities(im_phase, im_yfp)
    
    # update the all_intensities
    intensities_all = intensities_all + intensities
In [22]:
plt.hist(intensities_all,bins=20);

again, let's turn this into a function

In [23]:
def find_intensities_all(operator,repressor):
    # filename structure
    filename_phase = "./data/lacI_titration/" + operator + "_" + repressor + \
                    "_phase_pos_*.tif"
    filename_yfp = "./data/lacI_titration/" + operator + "_" + repressor + \
                    "_yfp_pos_*.tif"

    # grab all the images for this given strain
    phase_names = np.sort(glob.glob(filename_phase))
    yfp_names = np.sort(glob.glob(filename_yfp))

    # number of positions
    n_images = len(phase_names)

    # list to store intensity values for *all* the images
    intensities_all = []

    # loop through the images
    for i in range(n_images):
    
        # read in the image
        im_phase = skimage.io.imread(phase_names[i])
        im_yfp = skimage.io.imread(yfp_names[i])

        # get the intensities from from the YFP image
        intensities = find_intensities(im_phase, im_yfp)

        # update the all_intensities
        intensities_all = intensities_all + intensities
        
    return intensities_all

let's see if our function works

In [24]:
plt.hist(find_intensities_all("O3","R1740"))
Out[24]:
(array([136., 442.,  26.,   4.,   5.,   0.,   0.,   0.,   0.,   1.]),
 array([128.94642857, 145.45643688, 161.96644518, 178.47645349,
        194.98646179, 211.4964701 , 228.00647841, 244.51648671,
        261.02649502, 277.53650332, 294.04651163]),
 <a list of 10 Patch objects>)

Step 4: Get the fold-change for a given strain

In [25]:
def find_fold_change(operator, repressor):
    # mean intensities for "auto", "delta", and "Rxxx"
    mean_auto = np.mean(find_intensities_all(operator,"auto"))
    mean_delta = np.mean(find_intensities_all(operator,"delta"))
    mean_R = np.mean(find_intensities_all(operator,repressor))
    
    # calculate the fold change
    fold_change = (mean_R - mean_auto) / (mean_delta - mean_auto)
    
    return fold_change

test out our fold change function

In [26]:
find_fold_change("O1","R1220")
Out[26]:
0.005653994427639293