Assignment 9, Part 1: Analyze two-color yeast images

From Course Wiki
Revision as of 14:27, 13 November 2018 by Steven Wasserman (Talk | contribs)

Jump to: navigation, search
20.309: Biological Instrumentation and Measurement

ImageBar 774.jpg


Written questions

  1. If you didn't turn it in last time, use tfest or another method of your choosing to estimate parameters for a second-order transfer function model of the amplifier circuit you measured in the last assignment at various gain setting.
  2. What mathematical model did Mettetal, et. al. use for the yeast response network? Express the model in the following forms: transfer function, poles and zeros, single differential equation, and coupled differential equations. Express your answer in terms of the undamped natural frequency, $ \omega_0 $, damping ratio $ /zeta $, and/or damped natural frequency $ \omega_D $.
  3. Find an expression for the step response and plot it over a range of values of $ \omega_0 $ and $ /zeta $. A hand-drawn plot is fine, but you should probably look into MATLAB's step function.
  4. Plot the frequency response over a range of $ \omega_0 $ and $ /zeta $ values that includes over damped, critically damped, and under damped.

Yeast image analysis code

Start by downloading Assignment9TestData.mat. This is an example of the data you will be collecting in assignment 10 with an 8 minute valve oscillation period. The movie was collected using the same 20.309 microscope that you built in Assignment 8. During an oscillation experiment, the microscope first recorded an image with blue illumination followed by one with green illumination. It had a 40x objective and a 125 mm tube lens - remember that results in a 25X magnification.

The sample is similar to the one in Mettetal et al.: we are using S. Cerevisiae with the protein Hog1 fused to GFP (which is excited by blue light) and an mRNA binding protein (we'll call MCP) fused to tagRFP (which is excited by green light) to locate the nucleus. Our goal is to calculate the average Hog1 intensity in the nucleus (<GFP_nucleus>) vs. the cytosol (<GFP_cytosol>) as a function of the high/low valve state.

In MATLAB, display the contents of the structure in the data file:

twoColorYeastTest = 

  struct with fields:

                         Movie: [544×728×2×16 double]
                          Time: [16×2 double]
        ValveOscillationPeriod: 480
     BlueCameraGainAndExposure: [5 5000000]
    GreenCameraGainAndExposure: [15 5000000]
                    ValveState: [16×1 logical]

The data structure contains all the important information from the experiment:

  • A two-color movie, where the 3d dimension is the color (1 = Blue, 2 = Green), and the fourth is the frame number. Notice that twoColorYeastTest.Movie has already converted to a double, but it is not yet normalized by 4095.
  • The time stamp of each frame recorded, in seconds since the start of the experiment. Column 1 contains the timestamps of the Blue images, column 2 is for the green images.
  • The oscillation valve period, in seconds.
  • The blue and green image camera settings in the format [gain, exposure], where exposure is in microseconds.
  • The valve state (0 for low salt, 1 for high salt) at the time of each image. If you need to know the valve state more precisely, it can be calculated with the following conditional statement: $ sin(2 \pi t/T) \ge 0 $, where T is the valve oscillation period, and t is the time since the start of the experiment (this is how the state is determined in the software).


Pencil.png

Choose a single frame of the movie. Make a four panel figure and display: the blue image, the green image (both with a scale bar), and their corresponding image histograms. Include an appropriate title for each panel.


Analyze a frame

The first step in our analysis code is to write a function to extract relevant data from each movie frame. Our goal is to locate and identify individual cells and corresponding nuclei and then calculate the average Hog1 intensity in either region.

FindingCellsAndNucleiExample.png

Start with the following basic layout of the function:

function hog1Response = AnalyzeFrame(TwoColorFrame)

     % 1. estimate the background levels for each image
     % 2. find nuclei regions using green images
     % 3. find cell regions from blue images
     % 4. calculate the average GFP intensity (which is proportional to Hog1 concentration) inside 
     %    and outside the nucleus for each individual cell

end

As you work through this part of the Assignment, fill out the AnalyzeFrame function to incorporate each of the outlined steps.

Estimating background levels

BlueImageHistogram.png

It's tricky to get good images of proteins expressed by living yeast cells. Not only do we need to use very high exposure times (on the order of 5 s) to capture a dim signal, but we want to image cells over long times (minutes to hours). These two factors work against us to produce high background levels caused by camera noise (especially dark noise), light leakage from the room, or autofluorescence of the surrounding medium. Any drift in the background level over time can obscure our results. One way to get around this is to use the information contained within each image to estimate the background level.

In this particular movie, the yeast cells only occupy a small proportion of the image pixels. This can be seen in the image histogram, which contains a large low-intensity peak in both blue and green images. We will estimate the distribution of background pixels by fitting the histogram with a Gaussian function: $ p_1 \exp\left( -\frac{(x-p_2)^2}{2 p_3^2}\right) $.

The parameter $ p_2 $ represents the mean background level, while $ p_3 $ represents the width, or standard deviation of the background. Note that this method is based on the assumption that the intensity of the illumination is uniform (or, similarly that that the background distribution is the same everywhere in your image). This may not be the best assumption. Feel free to explore other ways of background correction in your analysis.

The following function can be used to estimate the mean and width of the pixel intensity distribution, assuming the background signal dominates:

function [backgroundMean, backgroundStd] = FitImageHistogram(frame)
    
    [counts, bins] = imhist(frame);
    figure
    plot(bins,counts)
    hold on
    
    gaussianFitFunction = @(p,x) p(1) * exp(-(x-p(2)).^2./(2 * p(3).^2));
    p1Guess = max(counts);
    p2Guess = bins(counts == p1Guess);
    [~, peakwidth] = min(abs(counts-p1Guess(1)/2));
    p3Guess = abs(bins(peakwidth) - p2Guess(1));
    beta0 = [p1Guess p2Guess p3Guess];
    
    fitParameters = nlinfit(bins,counts,gaussianFitFunction,beta0);
    
    plot(bins, gaussianFitFunction(fitParameters,bins));
    backgroundMean = fitParameters(2);
    backgroundStd = fitParameters(3);
    
    legend('image pixel data','best fit')
    
end


Pencil.png

Estimate the background mean and std for the blue and green color images of a single frame of the test movie.


Finding nuclei

We'll use a familiar method of locating the nuclei in our green image by setting a threshold to find the bright regions. In addition, we'll use a morphological opening operation (an erosion followed by a dilation) to remove hot pixels or other noise.

     nucleiMask = imbinarize(nucleiFrame,nucleiThreshold);
     nucleiMask = imopen(nucleiMask,strel('disk',2));

How should we choose the threshold? Since we just estimated the background level, one option is to use this information to pick a threshold well above the background level. Choosing a threshold greater than $ 3\sigma $ of the mean should result in eliminating 99% of the background pixels.

Finding cells

We could similarly use a threshold to identify cells in a blue image. The problem here, is that thresholding doesn't allow us to identify single cells within a clump of cells. Since we know that yeast cells are roughly spherical, we can use imfindcircles, which implements the Hough transform.

Here's a typical implementation:

     [cellCenters, cellRadii] = imfindcircles(cellFrame,[Rmin Rmax],'ObjectPolarity','bright','Sensitivity',0.95);

where Rmin and Rmax are the smallest and largest expected radii. How can you estimate Rmin and Rmax?


Loop through cells and calculate <GFP_cytosol> and <GFP_nucleus>

11/12/18 This section is still incomplete. stay tuned! -JS

A potentially useful helper function:

function Im = createMaskFromCircles(centers,radii,ImageSize)
    Im = zeros(ImageSize);
    [X, Y] = meshgrid(1:ImageSize(2), 1:1:ImageSize(1));

    for ii = 1:length(radii)
        isBright = (X-centers(ii,1)).^2+(Y-centers(ii,2)).^2 <= radii(ii,1)^2;
        isAlreadyFound = (Im == 1);
        Im = Im + (isBright & ~isAlreadyFound);
    end
    Im( Im > 1 ) = 1;
end

Eliminating bad cells