This tutorial is intended to demonstrate how to use MATLAB tools
to extend the functionality of
MicrobeTracker. A typical problem in the
analysis of microscope images is identification of large intracellular objects,
such as protein filaments, nucleoids, injected particles, or protein aggregates.
When big enough, these objects are very different from nearly
diffractionlimited spots and therefore require special treatment. The desired
propertis could be their number, size (area on an image), total intensity, or
centroid location. This tutorial demonstrates how to use standard MATLAB
functions perform (including some of the functions of the
Image Processing and
Mapping toolboxes) in
order to perform some of the mentioned tasks and to provide a framework for the
others.
The image set for this tutorial is of E. coli cells grown in M9 media
expressing HUCFP and imaged in phase contrast and fluorescence modes on an
agarose pad, so that their nucleoids are visible (courtesy of Dr. Manuel
Campos, JacobsWagner lab). The tutorials shows how to count the number of
nucleoids per cell, calculate their total and relative areas, as well as the
fluorescence intensity.
The data for the tutorials is not available by default. Please, make sure you
download the version of the MicrobeTracker Suite "with examples" from the
MicrobeTracker website.
The tutorial is divided into two sections: Pixelbased processing and
Subpixelresolution processing. The first one provides an easier mode of
analysis, sufficient in most cases. The second one uses some advanced Mapping
toolbox functions and is helpful when high precision is required or the cells
and intracellular objects are only a few pixels in size.
The tutorial assumes that the cell meshes are already obtained, they are
available in the meshes.mat file in the
microbetracker_nucleoid subfolder. This subfolder
also contains the phase contrast image used to generate the meshes
(phase.tif) and the fluorescence image used in this
tutorial to calculate the signal (fluo.tif). The
tutorial assumes that the user knows how to load the meshes into the MATLAB
workspace and how to run the included scripts. It is also recommended for the
user to be familiar with the basics of the Image Processing toolbox, which can
be achieved by going through the
Getting
Started tutorials on MATLAB's website. For reference, the images with
detected meshes are shown below (click on an image to see a zoomed version).
Pixelbased processing
1. Load the meshes and the fluorescence image into the MATLAB
workspace.
load('meshes.mat','cellList')
image = loadimagestack('fluo.tif');
Here MicrobeTracker's loadimagestack
function is used, capable of loading multiple images. To load a single image,
MATLAB's imread
function can be used instead.
2. Images in MATLAB can be stored either as arrays of unsigned
integers or as doubles. By default, a loaded image will be in the unsigned
integer format at the precision of the number matching that of the image file.
Many processing functions, however, require the image in double format, which
takes more space. When processing stacks of images, it is recommended to convert
them to doubles only one image at a time to save memory:
image = im2double(image);
3. Let's perform processing for the cell #2 in the dataset. The
dataset contains one frames, therefore the frame number will be 1. To reduce the
amount of text to write, get the variables for this particular cell we need: the
"box" identifying the area of the image containing the cell, and the mesh of the
cell (see Output Format section
for a description of these variables):
box = cellList{1}{2}.box;
mesh = cellList{1}{2}.mesh;
4. Crop the fluorescence image around the cell using
imcrop
function:
img1 = imcrop(image,box);
And convert the meshes to the cropped image coordinates:
mesh(:,[1 3]) = mesh(:,[1 3])box(1)+1;
mesh(:,[2 4]) = mesh(:,[2 4])box(2)+1;
5. Display the image and the mesh outline:
imshow(img1,[]);
hold on
plot(mesh(:,1),mesh(:,2),'y',mesh(:,3),mesh(:,4),'y')
set(gca,'position',[0 0 1 1])
Here the second argument of
imshow
function ([]) means automatic scaling of the image. In the second line,
hold on is
used to retain the image when the outline if drawn. The last line uses
set command
to change the axes
position and expand the image in the figure.
6. Convert the mesh to a polygon in order to use MATLAB's commands
operating with this format. The polygon consists of two variables  x and y
coordinates of the points on the outline of the cell:
x0 = [mesh(:,1);flipud(mesh(1:end1,3))];
y0 = [mesh(:,2);flipud(mesh(1:end1,4))];
7. Here is one example of using the polygon. Produce a mask of the cell using
poly2mask
command  a binary image with 1's corresponding to the pixels inside the cell
and 0's to the pixels outside:
cellmask = poly2mask(x0,y0,box(4)+1,box(3)+1);
8. The next step is to identify the nucleoid. One way to distinguish
the nucleoid (or another object) from the background is to find a threshold
intensity such that the regions of higher intensity would correspond to the
nucleoid, and the objects of lower intensity to the background. For this
MATLAB's graythresh
function can be used. To use this function, first normalize the image to span
the range of intensities from 0 to 1, since othervise it may fail:
img2 = img1min(img1(:));
img2 = img2/max(img2(:));
And then run graythresh function. Here only the pixels inside the cell are used
to calculate the threshold, they are obtained by using the cell mask as an index
for the image variable, i.e. img2(cellmask):
g = graythresh(img2(cellmask));
To get a mask corresponding to the nucleoid, get all the points with
intensities above the threshold:
nucleoidmask = img2>g;
There are alternatives to graythresh function to find objects in an image.
One popular alternative is applying a "Mexican hat" fitler to the image in
order to enhance the boundaries (use MATLAB's
fspecial
command with 'log' parameter to create the filter and
imfilter
to filter the image), and then take the areas below zero (sometimes below a
different threshold). This approach is beneficial when the object is assumed to
be of about uniform intensity and to have sharp boundaries, corrupted by
diffraction. In such case enhancing the boundaries is meaningful. In the case of
a nucleoid, where DNA density may decay very gradually, using boundaryenhancing
filters may result in rejection of weak areas neighboring brighter areas of the
nucleoid. In such cases, simple thresholding is prefered and therefore it was
chosen here.
9. Let's now visualize the nucleoid. For this use imshow command as
before. We'll combine the cell and the nucleoid masks to make the nucleoid white
on gray background of the cell and black background outside.
imshow(cellmask*0.5+nucleoidmask*0.5)
To display the cell outline, use the "hold on" command first, and then
visualize the outline polygon:
hold on
plot(x0,y0,'y')
10. Now let's cycle through all the cells in this example visualizing
the nucleoids in them. For this, open in the editor the script
get_nucleoid_pixilated_vis.m provided with this
example. I.e. type:
edit get_nucleoid_pixilated_vis
And then set a breakpoint in line 27 by clicking the dash after the line
number. After that, run the script by clicking Run button:
The script will display a cell and stop at the breakpoint. To continue to the
next cell, click Continue button or press F5.
11. Let's now calculate the number of nucleoids per cell. For this,
use MATLAB's regionprops
function to identify discontinuous regions in the nucleoid mask:
regstats = regionprops(nucleoidmask);
Then cycle through all regions picking those with the area larger than some
prefixed area areamin. On each step, if the region
is large enough, the counter nucleoidcount is
incremented:
areamin = 50;
nucleoidcount = 0;
for i=1:length(regstats)
if regstats(i).Area>=areamin
nucleoidcount = nucleoidcount+1;
end
end
nucleoidcount
A script that cycles through each cell, calculates the number of nucleoids,
and then displays a histogran of the number of nucleoids per cell is provided
with this example and is called
get_nucleoid_count_pixilated.m. Now run it:
get_nucleoid_count_pixilated
12. Let's now calculate the nucleoid areas. With a small modification
of the above script, we now append the areas of the nucleoids we were counting
before:
nucleoidareaarray = [];
areamin = 50;
for i=1:length(regstats)
if regstats(i).Area>=areamin
nucleoidareaarray = [nucleoidareaarray regstats(i).Area];
end
end
Here is another script provided with the example that displays a histogram of
nucleoid areas:
get_nucleoid_area_pixilated
13. Now, let's calculate the fraction of the cell occupied by
nucleoids. This time we don't need the regoinprops function. We calculate the
areas of both the cell and the nucleoid as sums of pixel values for the masks,
as the masks have value 0 for the pixels outside and 1 for the pixels inside
the object:
cellarea = sum(cellmask(:));
nucleoidarea = sum(nucleoidmask(:));
To get a fraction, simply divide one by the other:
nucleoidarea/cellarea
Another script provided with the example displays a histogram of the relative
areas occupied by nucleoids:
get_nucleoid_relative_area_pixilated
14. Finally, let's calculate the average intensity of the signal
inside cells and inside nucleoids for every cell, and display as two histograms.
For this, an extra step of removing the background from the fluorescence image
is required. Run MicrobeTracker, load phase contrast and fluorescence images,
load parameters, select All radio button on the
Background subtraction
panel and click Subtract bgrnd, then click Save vis. signal and
save the buckgroundsubtreacted fluorescence image to a new file (here
fluo_subtracted.tif was used). See the image below
for a sequence of operations in MicrobeTracker:
To calculate the total signal in the nucleoid sum the image pixels in the
mask:
sum(img1(nucleoidmask))
And divide them by the total area, i.e. for the nucleoids and for the cells:
nucleoidint = sum(img1(nucleoidmask))/sum(nucleoidmask(:));
cellint = sum(img1(cellmask))/sum(cellmask(:));
The rest of the procedure is similar to the procedures described above. The
whole script calculating the intensities (total signal divided by the area) in
nucleoids and in the cells (including nucleoids) is called
get_nucleoid_intensities_pixilated.m. Now execute
it:
get_nucleoid_intensities_pixilated
Subpixelresolution processing
1. The beginning of subpixelresolution processing is the same up to
the point of determining the threshold. I.e. for cell #1 on frame #1:
load('meshes.mat','cellList')
image = im2double(loadimagestack('fluo.tif'));
box = cellList{1}{1}.box;
mesh = cellList{1}{1}.mesh;
x0 = [mesh(:,1);flipud(mesh(1:end1,3))]box(1)+1;
y0 = [mesh(:,2);flipud(mesh(1:end1,4))]box(2)+1;
img1 = imcrop(image,box);
mask = poly2mask(x0,y0,box(4)+1,box(3)+1);
img2 = img1min(img1(:));
img2 = img2/max(img2(:));
g = graythresh(img2(mask));
2. Then we find the outline of the nucleois or nucleoids in the cell.
First, set the value of the pixels outside of the cell to zero so that no
objects are detected outside. However, a fraction of a pixel may be inside of
the cell even though its center is outside. Therefore, we first expand the cell
mask by one pixel outwards:
expandedmask = imdilate(mask,strel('square',3));
And then set the values of the outside pixels to zero (the ~ sign means
negation):
img2(~expandedmask)=0;
Now we can determine the contours at subpixel resolution by using MATLAB's
contourf
function. This function uses linear interpolation between pixels to get the
outline as precise as the information in the image allows, not limited to the
pixel size:
c = contourf(img2,[g g]);
The format of the output will be explained below.
2. Now lets parse the contour structure obtained with the contourf
command. The structure consists of the lines of x and y coordinates of
individual contours, separated by the values of the polygons (equal to g here)
and the number of vertices in each polygon. See
contour
matrix description for details. To parse the structure, we can go through,
extracting the individual polygons (with coordinates of the vertices x and y) on
the way:
ind = 1;
while ind<size(c,2)
ctr = c(:,ind+1:ind+c(2,ind))';
ind = ind+c(2,ind)+1;
x1 = ctr(:,1);
y1 = ctr(:,2);
end
3. In order to take only the portions of the polygons inclosed in the
cell outline, we can use MATLAB's
polybool
function:
[x2,y2] = polybool('intersection',x0,y0,x1,y1);
Note that polybool function assumes that the contours have clockwise
orientation, and though it works, it displays a warning otherwise. To convert
counterclockwiseoriented contours to clockwise orientation, use MATLAB's
poly2cw
function (this step is not shown here, but included in the final script).
Let's display the (x1,y1) and (x2,y2) polygons to see what we achieved so
far. Copy and paste the code below to MATLAB's command window (assuming the
contour cstructure has been already generated). The cell outline will be
plotted in yellow color, the original polygons in green, and their portions
inside the cell in red:
figure
imshow(img1,[])
set(gca,'pos',[0 0 1 1])
hold on
plot(x0,y0,'y')
ind = 1;
while ind<size(c,2)
ctr = c(:,ind+1:ind+c(2,ind))';
ind = ind+c(2,ind)+1;
x1 = ctr(:,1);
y1 = ctr(:,2);
[x2,y2] = polybool('intersection',x0,y0,x1,y1);
plot(x1,y1,'g')
plot(x2,y2,'r')
end
4. Now separate the outlines of the nucleoids from the outlines of the
included empty areas. This can be done by determining the orientation of the
polygons included in the contour cstructure: the outer polygons will have
clockwise orientation, while the internal ones  counterclockwise. To determine
the orientation or a polygon, MATLAB's
ispolycw
function can be used (or a faster MicrobeTracker's isContourClockwise function).
Run included get_nucleoid_subpixel_vis.m script
in debug mode to visualize the nucleoids in every cell onebyone. For this,
open the script in the editor, put a breakpoint in line 38, and run the script.
The nucleoids will be shown as cyan patches, the "holes" in them will be colored
yellow.
5. Finally, let's calculate the relative areas occupied by nucleoids.
For this, go through the nucleoid regions and add the areas of those with
clockwise orientation (external), substracting those with counterclockwise
orientation (internal):
cellarea = polyarea(x0,y0);
nucleoidarea = 0;
ind = 1;
while ind<size(c,2)
ctr = c(:,ind+1:ind+c(2,ind))';
ind = ind+c(2,ind)+1;
x1 = ctr(:,1);
y1 = ctr(:,2);
[x2,y2] = polybool('intersection',x0,y0,x1,y1);
if ispolycw(x1,y1)
nucleoidarea = nucleoidarea+polyarea(x2,y2);
else
nucleoidarea = nucleoidareapolyarea(x2,y2);
end
end
The relative area will be a ratio of nucleoidarea
to cellarea:
nucleoidarea/cellarea
The whole script calculating this ratio for every cell and displaying it in
for of a hystogram is called
get_nucleoid_relative_area_subpixel.m. Now execute
it:
get_nucleoid_relative_area_subpixel
It is not difficult to get individual nucleoid areas. However, the problem of
precisely integrating the signal inside arbitrary polygons is significantly more
challenging. Those who are interested can refer to
getOneSignalC function in MicrobeTracker's code
itself.
