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 0x1c1bde0a90>
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 0x1c1cbad9b0>
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 0x1c1d4fa208>

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 0x1c1deed2b0>

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 0x1c1df97780>
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 0x1c1e3b9470>

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 0x1c1e91a0b8>

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 0x1c1f9d77f0>

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

Step 5: Analyze the whole dataset

In [27]:
# names of the strains we are dealing with
operators = ["O1","O2","O3"]
repressors = ["R22","R60","R124","R260","R1220","R1740"]
repressor_counts = [22,60,124,260,1220,1740]

# number of operators and repressor
n_operators = len(operators)
n_repressors = len(repressors)

# array to store our fold change values in 
fc_values = np.zeros([n_operators, n_repressors])

# first, loop through operators
for i in range(n_operators):
    
    # second, loop through number of repressors
    for j in range(n_repressors):
        
        # grab the operator / repressor pair
        operator = operators[i]
        repressor = repressors[j]
        
        # print output to see progress through for loop
        print(operator, repressor)
        
        # find the find-change
        fc_values[i,j] = find_fold_change(operator, repressor)

O1 R22
O1 R60
O1 R124
O1 R260
O1 R1220
O1 R1740
O2 R22
O2 R60
O2 R124
O2 R260
O2 R1220
O2 R1740
O3 R22
O3 R60
O3 R124
O3 R260
O3 R1220
O3 R1740

Step 6: Plot the theory curves

Plotting the three states over increasing repressors

In [28]:
# parameters 
dE_R = -9 # binding energy of repressor to DNA, units of kBT
dE_P = -6 # binding energy of polymerase to DNA, units of kBT
P = 1000 # number of RNAP per cell
N_NS = 5E6

# range of repressor copy numbers to look at
R = np.logspace(0,4,100)

# weights
empty_weight = np.ones(len(R))
R_bound_weight = (R / N_NS) * np.exp(-dE_R)
P_bound_weight = np.ones(len(R)) * (P / N_NS) * np.exp(-dE_P)
tot_weight = empty_weight + R_bound_weight + P_bound_weight 

# compute the probabilities
empty_prob = empty_weight/tot_weight
R_bound_prob = R_bound_weight/tot_weight
P_bound_prob = P_bound_weight/tot_weight
tot_prob = empty_prob + R_bound_prob + P_bound_prob

# plot
plt.plot(R, empty_prob)
plt.plot(R, R_bound_prob)
plt.plot(R, P_bound_prob)
plt.plot(R, tot_prob, "k")

# plot labeling
plt.xscale("log")
plt.xlabel("repressor copy number")
plt.ylabel("probability")
plt.legend(["empty", "repressor bound", "polymerase bound", "total"])
Out[28]:
<matplotlib.legend.Legend at 0x1c1fdf3e48>

plot the theory curves for fold change

In [29]:
# repressor binding energies in kBT units
dE = [-15.3, -13.9, -9.3]

# colors for plotting
colors = ["k", "g", "b"]

# calculate and plot the theoritcal prediction
for i in range(len(dE)):
    fold_change = 1 / (1 + R / N_NS * np.exp(-dE[i]))
    plt.loglog(R, fold_change, color=colors[i])
    plt.loglog(repressor_counts, fc_values[i,:], '.', color=colors[i])
    
plt.xlabel("repressor copy number")
plt.ylabel("fold-change")
plt.legend(["O1 theory", "O1 data", \
            "O2 theory", "O2 data", \
            "O3 theory", "O3 data"])
Out[29]:
<matplotlib.legend.Legend at 0x108d37cf8>

Step 7: Extra credit version -- now including induction

In [30]:
# energies, in kBT
dEI = -1 # energy of ligand bound to inactive
dEA = -5 # energy of ligand bound to active 
dES = 3 # energy of the active state (regardless of ligand)

# ligand concentration range to look at
ligand_conc = np.logspace(-3, 2, 100)

# weights
# I for inactive, A for active, U for unbound, B for bound 
IU_w = np.ones(len(ligand_conc))
IB_w = ligand_conc * np.exp(-dEI)
AU_w = np.exp(-dES) * np.ones(len(ligand_conc))
AB_w = ligand_conc * np.exp(-dES) * np.exp(-dEA)
tot_w = IU_w + IB_w + AU_w + AB_w

# convert the weights to probabilities
IU_p = IU_w / tot_w
IB_p = IB_w / tot_w
AU_p = AU_w / tot_w
AB_p = AB_w / tot_w

# plot
plt.plot(ligand_conc, IU_p, "r--")
plt.plot(ligand_conc, IB_p, "r-")
plt.plot(ligand_conc, AU_p, "k--")
plt.plot(ligand_conc, AB_p, "k-")

plt.xscale("log")
plt.legend(["inactive unbound", "inactive bound", "active unbound", "active bound"])
Out[30]:
<matplotlib.legend.Legend at 0x1c206d8c50>