First move to the directory where the python files for image analysis are stored.

In [1]:
pushd "../../../Microscope_image_processing_code/microscopyimageprocessing"
C:\Users\Nicholas Sherer\Box\N_Sherer_Thesis_Materials\Microscopy_Experiments\Microscope_image_processing_code\microscopyimageprocessing
Out[1]:
['~\\Box\\N_Sherer_Thesis_Materials\\Microscopy_Experiments\\Microscope_Images\\30 November 2018\\NS001 2000 IPTG 6 aTc TIFFS']

Then import the necessary python libraries (including our custom collections of helper functions) and set up our matplotlib plotting options.

In [2]:
import copy
from datetime import datetime
import os

import numpy as np
import matplotlib as mpl
mpl.rcParams["axes.grid"] = False
mpl.rcParams["image.cmap"] = "viridis"
mpl.rcParams["image.interpolation"]="nearest"
import matplotlib.pyplot as plt
%matplotlib inline

import IPython.display as ipyd

import skimage.morphology as skmo
import skimage.measure as skme
import skimage.io as skio
import skimage.transform as sktr

import illuminationinterpolation as illint
import segmentation as mseg
import visualization as mvis
import savenbs

NOTE: below command allows us to go back to the directory in which the notebook is stored in. Simply drop the notebook along with the data in the same folder and you're good to go.

In [3]:
popd
C:\Users\Nicholas Sherer\Box\N_Sherer_Thesis_Materials\Microscopy_Experiments\Microscope_Images\30 November 2018\NS001 2000 IPTG 6 aTc TIFFS
popd -> ~\Box\N_Sherer_Thesis_Materials\Microscopy_Experiments\Microscope_Images\30 November 2018\NS001 2000 IPTG 6 aTc TIFFS

Loading microscope images are located and naming the save folder

First we need to give the filenames and paths to where our images of beads and e. coli are located so they can be loaded. Then we need to give glob patterns so we can sort the various image channels from each other. Finally, we give a name to the folder we'll be saving results in. We append the current time (in UTC) to the save folder name to avoid accidentally overwriting past analysis.

In [4]:
beads_folder_name = ""
beads_name = "beads"
beads_pc = skio.ImageCollection(beads_folder_name + beads_name + "*c1*c2.tif")
beads_TIRF_561 = skio.ImageCollection(beads_folder_name + beads_name + "*c2*c2.tif")
beads_BF = skio.ImageCollection(beads_folder_name + beads_name + "*c2*c1.tif")

cell_folder_name =""
cell_name = "ns001"
ecoli_TIRF_561 = skio.ImageCollection(cell_folder_name + cell_name + "*c2*c2.tif")
ecoli_pc = skio.ImageCollection(cell_folder_name + cell_name + "*c1*c2.tif")

save_folder = "standard analysis " + repr(datetime.utcnow())[9:]
os.mkdir(save_folder)

Aligning the DS-Fi2 camera to the DU-897 camera

First we need to use the images of beads for aligning the camera and finding the transformations relating the coordinates between different types of images.

We first inspect the bead images to make sure they've loaded and they look good.

In [5]:
mvis.inspectImages([beads_BF, beads_TIRF_561, beads_pc], (16,10))

Then find the background illumination in the brightfield and phase contrast channels and normalize them by dividing each image by the background illumination.

In [6]:
beads_BF_bg = mseg.findMedianBg(beads_BF)
beads_pc_bg = mseg.findMedianBg(beads_pc)
beads_BF_n = [image / beads_BF_bg for image in beads_BF]
beads_pc_n = [image / beads_pc_bg for image in beads_pc]
In [7]:
plt.plot([np.mean(image) for image in beads_BF],'*')
plt.xlabel('FOV')
plt.ylabel('mean brightness')
Out[7]:
Text(0, 0.5, 'mean brightness')

We see not much pattern in mean image brightnesses.

Next before the phase contrast and TIRF images can be aligned, we need to find the beads. This function does it for brightfield, but a simple threshold will work for phase contrast or TIRF. Always inspect the masks to make sure they look good. It's ok if a few beads are missed.

In [8]:
threshes = [.9*np.mean(image) for image in beads_BF_n]
bead_BF_masks = [mseg.findBeadsBF(image,thr) for image, thr in zip(beads_BF_n,threshes)]
beads_BF_masks = mseg.removeNonCirles(bead_BF_masks,10)
mvis.inspectImages([beads_BF_n, beads_BF_masks],figsize=(14,7))

The phase contrast channel is a lot easier to find the masks for. A high threshold usually works. Inspect these masks next to the brightfield masks. We'll need a good field of view for alignment.

In [9]:
bead_pc_masks = [image>2.0*np.mean(image) for image in beads_pc_n]
bpcm = [skmo.remove_small_holes(mask,min_size=3000) for mask in bead_pc_masks]
bpcm = [skmo.remove_small_objects(mask, min_size=1000) for mask in bpcm]
mvis.inspectImages([beads_pc_n, bpcm, beads_BF_masks],(14,7))
E:\Anaconda3\envs\imageprocessing_windows\lib\site-packages\skimage\morphology\misc.py:206: UserWarning: the min_size argument is deprecated and will be removed in 0.16. Use area_threshold instead.
  warn("the min_size argument is deprecated and will be removed in " +

Once you've picked a field of view, the findRegionCenters helper function will find the centers of the beads in a mask.

In [10]:
index2=1
BF2pc_align_mask = bead_BF_masks[index2]
PC_align_mask = bpcm[index2]
BF2pc_centers = mseg.findRegionCenters(BF2pc_align_mask)
PC_centers = mseg.findRegionCenters(PC_align_mask)

You'll need to loop through the next two cells until the lines drawn between the two images are matching beads up correctly. The first cell is used to add, removem and rearrange the order of bead centers from the lists of centers to be matched between cameras.

In [11]:
BF2pc_centers_cl = BF2pc_centers[[0,3,4,5,7,8,9]]
PC_centers_cl = PC_centers
In [12]:
mvis.showKeypointpairs(BF2pc_align_mask, PC_align_mask, BF2pc_centers_cl, PC_centers_cl)

Once the beads are paired up well, skimage can easily find the geometric transform between them.

It's good to print out the transform because we may need to adjust it by hand later down. We also make a copy of it so we don't adjust the original transform skimage found

In [13]:
camera_transform_p = sktr.SimilarityTransform()
camera_transform_p.estimate(BF2pc_centers_cl,PC_centers_cl)
beads_PC_warped = mseg.warpIm2Im(beads_pc_n[index2],beads_BF_n[index2],camera_transform_p)
mvis.showImages([beads_PC_warped, beads_BF_n[index2]],(14,7))
total_transform = sktr.SimilarityTransform()
total_transform.params = copy.deepcopy(camera_transform_p.params)
print(total_transform.params)
[[ 2.36803804e+00 -1.63752257e-01  1.65259837e+01]
 [ 1.63752257e-01  2.36803804e+00 -2.93864095e+02]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]

Finding the illumination field in 561

Now that we've found the alignment between the different types of images; the beads can be used as a reference to find the illumination field which can correct intensities for the unevenness of the TIRF illumination field.

First we find masks for each bead. We can use brightfield masks for this or warp the phase contrast masks into brightfield coordinates. Using the phase contrast masks will lose some information on the illumination field but we only consider e. coli we see in phase contrast anyways so that lost information is irrelevant. Making masks from the TIRF directly is difficult and error prone.

In [14]:
def background(mask, radius):
    return np.logical_not(skmo.binary_dilation(mask,selem=skmo.disk(radius)))

background_masks = [background(mask, 60) for mask in beads_BF_masks]

backgrounds_masked_561 = np.ma.masked_array(np.array(beads_TIRF_561), mask=np.array(np.logical_not(background_masks)))
median_561_background = np.ma.median(backgrounds_masked_561, axis=0)
if save_folder:
    with savenbs.cd(save_folder):
        np.savez_compressed('median_561_background', median_561_background.data, allow_pickle=False)
plt.imshow(median_561_background)
Out[14]:
<matplotlib.image.AxesImage at 0x2aeb92245f8>
In [15]:
beads_TIRF_561_bg_sub = [image - median_561_background for image in beads_TIRF_561]
beads_label = [skme.label(mask) for mask in beads_BF_masks]
beads_rprops_561 = [skme.regionprops(label, TIRF) for label, TIRF in zip(beads_label, beads_TIRF_561_bg_sub)]
beads_plist_561 = mseg.properties2list(beads_rprops_561,['mean_intensity','area','centroid'])

print('There are', len(beads_plist_561['mean_intensity']), 'beads.')
print('The mean 561 intensity of the beads is', '%.0f' % np.mean(beads_plist_561['mean_intensity']), '.')
print('The standard deviation of 561 intensity of the beads is', 
      '%.0f' % np.std(beads_plist_561['mean_intensity']), '.')
print('The mean area of the beads is', '%.0f' % np.mean(beads_plist_561['area']), '.')
print('The standard deviation of area of the beads is', '%.0f' % np.std(beads_plist_561['area']), '.')
There are 212 beads.
The mean 561 intensity of the beads is 20375 .
The standard deviation of 561 intensity of the beads is 8744 .
The mean area of the beads is 506 .
The standard deviation of area of the beads is 20 .
In [16]:
bead_gauss_dfunc_561 = illint.createGaussianDistFunc(beads_plist_561['centroid'],
                                                     beads_plist_561['mean_intensity'],40)
bead_TIRF_bg_gaus_561 = illint.evaluateDistsBox(512,512,bead_gauss_dfunc_561)
bead_TIRF_bg_gaus_561 = bead_TIRF_bg_gaus_561 / np.mean(bead_TIRF_bg_gaus_561)
if save_folder:
    with savenbs.cd(save_folder):
        np.savez_compressed('inferred_561_illumination', bead_TIRF_bg_gaus_561, allow_pickle=False)
In [17]:
plt.imshow(bead_TIRF_bg_gaus_561)
bg_max_561 = np.max(bead_TIRF_bg_gaus_561)
bg_min_561 = np.min(bead_TIRF_bg_gaus_561)
print('The brightest relative intensity is', '%.2f' % bg_max_561, '.')
print('The dimmest relative intensity is', '%.2f' % bg_min_561, '.')
print('This is a factor of', '%.2f' % (bg_max_561/bg_min_561), 'between the brightest and dimmest areas.')
The brightest relative intensity is 1.79 .
The dimmest relative intensity is 0.30 .
This is a factor of 5.86 between the brightest and dimmest areas.
In [18]:
beads_TIRF_561_n = [image/bead_TIRF_bg_gaus_561 for image in beads_TIRF_561_bg_sub]
beads_rprops_561_n = [skme.regionprops(label, TIRF) for label, TIRF in zip(beads_label, beads_TIRF_561_n)]
beads_plist_561_n = mseg.properties2list(beads_rprops_561_n,['mean_intensity','area','centroid'])
if save_folder:
    with savenbs.cd(save_folder):
        np.save('beads_561_intensities', beads_plist_561_n['mean_intensity'], allow_pickle=False)
n_beads = len(beads_plist_561_n['mean_intensity'])
plt.hist(beads_plist_561_n['mean_intensity'],bins=np.sqrt(n_beads).astype('int32'));
bead_mean_561_n = np.mean(beads_plist_561_n['mean_intensity'])
bead_stderr_561_n = np.std(beads_plist_561_n['mean_intensity'])/np.sqrt(n_beads)
print('After illumination correction, the mean intensity of the beads is',
      '%.0f' % bead_mean_561_n,'+-', '%.0f' % bead_stderr_561_n)
After illumination correction, the mean intensity of the beads is 20362 +- 161

Identifying the E. coli in phase contrast

We use the phase contrast images of E. coli to separate the cells from the background. We also use the illumination field gained from the beads to correct the illumination in the E. coli TIRF images.

First, we inspect the imported images of E. coli.

In [19]:
mvis.inspectImages([ecoli_pc, ecoli_TIRF_561],(18,8))

This next function does the illumination correction for the e. coli phase contrast images and also nonlocal means denoising (using opencv library) on the phase contrast. This denoising is important to getting a good result from thresholding to find the e. coli. The denoising will generally give a warning about loss of precision upon converting float64 to uint8. Do not worry about this, phase contrast intensities are not used for quantification so it doesn't matter.

Once we've normalized the images, it's time to threshold the phase contrast images. The automated thresholding followed by cleaning up any objects too small to be e. coli works pretty well, but goes wonky occasinally so we'll have to manually drop any images with bad segmentation.

In [20]:
ecoli_pc_norm = mseg.normAndDenoisePc(ecoli_pc)
b_masks_cl = [mseg.thresholdMask(image,min_size=300) for image in ecoli_pc_norm]
if save_folder:
    with savenbs.cd(save_folder):
        np.savez_compressed('ecoli_pc_norm', ecoli_pc_norm, allow_pickle=False)

Of course, we inspect the masks as we inspect everything else. The index on the slider of masks that look bad should go into the indices to blank function of the cell after.

In [21]:
mvis.inspectImages([b_masks_cl,ecoli_pc_norm],(24,16))

Then we drop any images with bad segmentation.

In [22]:
indices_to_blank = []
for index in indices_to_blank:
    b_masks_cl[index] = np.zeros_like(b_masks_cl[index])
if save_folder:
    with savenbs.cd(save_folder):
        np.savez_compressed('phase_contrast_masks_cleaned', b_masks_cl, allow_pickle=False)

Inspecting aligment and making small corrections to misalignment.

Next I use the masks for e. coli from phase contrast and transform them to align to the TIRF and see how well they line up. I then tweak the parameters of the transformation by hand using the sliders if the alignment is unsatisfactory to my eyes.

In [23]:
camera_transform_ec = copy.deepcopy(total_transform)
tabwidget, interaction = mvis.adjustAlignment(ecoli_TIRF_561, b_masks_cl, camera_transform_ec)
ipyd.display(tabwidget, interaction)

The transform looks good.

In [24]:
if save_folder:
    with savenbs.cd(save_folder):
        np.save('camera_to_camera_transform', camera_transform_ec.params, allow_pickle=False)
print(camera_transform_ec.params)
[[ 2.36803804e+00 -1.63752257e-01  1.65259837e+01]
 [ 1.63752257e-01  2.36803804e+00 -2.93864095e+02]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]

Next we make labeled arrays from the phase contrast masks and then only after labeling, warp them into the TIRF coordinates (otherwise we accidentally might fuse separated cells since there are more pixels per area in phase contrast than in TIRF).

In [25]:
ec_labels = [skme.label(mask) for mask in b_masks_cl]
ec_labels_w = [mseg.warpIm2Im(label,ecoli_TIRF_561[0],camera_transform_ec) for label in ec_labels]
ec_pc_w = [mseg.warpIm2Im(image,ecoli_TIRF_561[0],camera_transform_ec) for image in ecoli_pc_norm]

Now that we've labeled the e. coli we can subtract the background and normalize the TIRF images. I hit some fields of view below that were too crowded to find a background near some e. coli. I'm going to run it again, but catch the runtime error and drop any fields of view where this is true. Hopefully that still leaves enough fields of view to be workable.

In [26]:
ecoli_TIRF_561_norm = []
ec_labels_w_success = []
failures = []
i=0
for image, label in zip(ecoli_TIRF_561, ec_labels_w):
    try:
        corrected = mseg.subtract_pad_bg(image, label, 5, 10)/bead_TIRF_bg_gaus_561
        ecoli_TIRF_561_norm.append(corrected)
        ec_labels_w_success.append(label)
    except RuntimeError:
        ecoli_TIRF_561_norm.append(np.zeros_like(image))
        ec_labels_w_success.append(np.zeros_like(label))
        failures.append(i)
    i = i + 1
if save_folder:
    with savenbs.cd(save_folder):
        np.savez_compressed('ecoli_TIRF_561_norm', ecoli_TIRF_561_norm, allow_pickle=False)
        np.savez_compressed('ecoli_labels_success', ec_labels_w_success, allow_pickle=False)
print(failures)
[]
In [27]:
mvis.inspectImages([ec_labels_w_success, ecoli_TIRF_561_norm],(24,16))

Determining the mean 561 signal (brightness per pixel in a cell)

Now the only thing that remains is to use the warped phase contrast masks and normalized TIRF images to find the distribution of e. coli intensities.

In [28]:
ec_rprops_561 = [skme.regionprops(label,TIRF) for label, TIRF in zip(ec_labels_w_success, ecoli_TIRF_561_norm)]
ec_intensities_561 = mseg.properties2list(ec_rprops_561, ['mean_intensity'])['mean_intensity']
n_ec = len(ec_intensities_561)

plt.figure()
urg = plt.subplot(111)
urg.hist(ec_intensities_561,bins=np.sqrt(n_ec).astype('int32'));

ec_meanintensity_561 = np.mean(ec_intensities_561)
ec_stderrintensity_561 = np.std(ec_intensities_561)/np.sqrt(n_ec)
print('The mean intensity of the e. coli is', '%.0f' % ec_meanintensity_561,'+-', '%.0f' % ec_stderrintensity_561)
The mean intensity of the e. coli is 325 +- 9

Comparing the signal to an objective standard

To make the signal meaningful for comparison, we have to compare it to the brightness of the beads on the same slide.

In [29]:
ec_intensities_n_561 = ec_intensities_561 / bead_mean_561_n
if save_folder:
    with savenbs.cd(save_folder):
        np.save('ecoli_561_intensities_normalized', ec_intensities_n_561, allow_pickle=False)

plt.hist(ec_intensities_n_561,bins=np.sqrt(n_ec).astype('int32'));

ec_meanintensity_n_561 = ec_meanintensity_561 / bead_mean_561_n
ec_meanintensity_n_lefterr_561 = (ec_meanintensity_n_561)/(bead_mean_561_n + bead_stderr_561_n)
ec_meanintensity_n_righterr_561 = (ec_meanintensity_n_561)/(bead_mean_561_n - bead_stderr_561_n)
print('The ecoli mean intensity is estimated to be', '%.5f,' %(ec_meanintensity_n_561),
      'with left and right standard erors of ', '(%.5f,' %(ec_meanintensity_n_lefterr_561), 
      '%.3f).' %(ec_meanintensity_n_righterr_561))

ec_stdintensity_n_561 = np.std(ec_intensities_n_561)
print('The standard deviation over the mean is', '%.3f.' %(ec_stdintensity_n_561/ec_meanintensity_n_561))
The ecoli mean intensity is estimated to be 0.01597, with left and right standard erors of  (0.00000, 0.000).
The standard deviation over the mean is 0.696.
In [ ]: