2D SDS Electrophoresis

The picture is from an online publication at this link: http://www.fao.org/docrep/008/y5970e/y5970e07.htm

Gel_Detect – detect and cluster spots on an 2D SDS gel

This demo shows how the image processing toolbox together with matrix operations and statistics toolbox can be used to process an SDS page gel. SDS -> sodium dodecyl sulfate polyacrylamid gel. SDS is an anionic detergent that denatures secondary structure of proteins. It is incubated with the protein sample so that the propagation of each protein just dependends on its molecular mass. In 2D electrophoresis, proteins are separated here in two dimensions, once by their isoelectric point and once by their mass/charge ratio.

The challenge in analysing the gel images is in the first instance to find out which spots are interesting resp. important and secondly to map spots to proteins. Sometimes this is just done with an overlay, to compare two samples or one sample with a known reference mix. This demo shows how we can determine the spots automatically and build clusters of them. The mapping has been left out – but could easily be done in addition.

Contents

• Read the image, convert it to grayscale
• Let the user crop the image
• Create a spot mask
• Label the spots in the image.
• Count the number of spots in the Gel
• Get region properties: centroid, bounding box and area
• Fetch centroids (mass centers) and bounding boxes of the spots
• Work with spot centroids, sizes and intensities
• Prepare spot data for clustering
• Get and display spot clusters via kmeans
• Overlay spot clusters with original image

Read the image, convert it to grayscale

gelGrayOrig = rgb2gray(imread('gel.tif'));
figure, imshow(gelGrayOrig), title('original image');

Let the user crop the image

Here the user can define the borders of the image. This is important for later measurements of the spot distances in pH gradient and SDS migration (size of the molecule)

gelGray = imcrop;
figure, imshow(gelGray), title('cropped image');

Create a spot mask

Pixels with an intensity value above a certain threshold can be set to zero. With this mechanism we can determine a mask that contains only the spots on the gel. First, copy the image, then zero out all values that are too “white”.

The sensitivity to spots depends heavily on the choosen threshold and also on the contrast of the image. It might be necessary to sharpen the picture and/or let the user choose the threshold manually and interactively. This can be achieved by wrapping the whole program in a GUI (easily possible with ML GUI building facilities), creating a slider for the threshold which has a callback that overlays the spot mask with the image for the actually chosen threshold.

gelMask = gelGray;
wTheta = 160;
gelMask(gelGray > wTheta) = 0;
gelIntens = gelMask;                % keep the intensities.
                                    % We want to refer to them later
gelMask(gelGray <= wTheta) = 255;
gelMask = im2bw(gelMask);           % convert to binary image
figure, imshow(gelMask), title('spot mask (gelMask)');

Label the spots in the image.

This is helpful for detection of different spots. Each connected area in the image (area of adjacent one’s) is labeled with a different number. The number of distinct labels equals the number of putative spots.

gelMaskLabel = bwlabel(gelMask);
figure, imshow(label2rgb(gelMaskLabel)), title('labeled gel mask (gelMaskLabel)');

Count the number of spots in the Gel

nSpots = max(max(gelMaskLabel));
disp(['Total number of spots: ' num2str(nSpots)]);

Total number of spots: 232

Get region properties: centroid, bounding box and area

The regionprops function returns several properties of labeled images, which can be defined as named arguments.

gelRegionProps = regionprops(gelMaskLabel, 'Area', 'BoundingBox', ...
    'Centroid', 'MajorAxisLength', 'MinorAxisLength');

Fetch centroids (mass centers) and bounding boxes of the spots

The bounding box contains the elements upper left corner of the spot (x y coordinates) and width (x y)

gelSpotCentroids = [gelRegionProps.Centroid];
gelSpotBBox = [gelRegionProps.BoundingBox];

Work with spot centroids, sizes and intensities

At this point we have a lot of spots with eventually different sizes, intensities and positions. The task is now to find the “interesting” spots – e.g. the ones that represent the protein we are looking for or the ones that show important facts about the sample. Eventually, this can be achieved with pattern recognition approaches. As important factors we would consider the bounding box (as a measure of spot size), the major/minor axis length (as a measure of spot extension), the area, the position and the intensity. Good factors would be ratios of area/bbox, area/(major,minor axis length).

Prepare spot data for clustering

Goal is a matrix that describes several spot properties. This can be used for some further analysis.

% compute bounding box diagonals
spotBBdiag = zeros(nSpots, 1);
j = 1;
for i = 3:4:length(gelSpotBBox)
    spotBBdiag(j) = sqrt(gelSpotBBox(i)^2 + gelSpotBBox(i+1)^2);
    j = j + 1;
end% get average intensities for each spot
spotAvgIntens = zeros(nSpots, 1);
for i = 1:nSpots
    spotAvgIntens(i) = mean(gelGray(find(gelMaskLabel == i)));
end
clear i j;

% create spot properties matrix. The rows are:
% Area Intensity MajorAxis MinorAxis BoundingBoxDiagonal
spotProperties = [[gelRegionProps.Area]' spotAvgIntens ...
    [gelRegionProps.MajorAxisLength]' ...
    [gelRegionProps.MinorAxisLength]' spotBBdiag];
spotPropNames = {'Area' 'Average Intensity' 'Major axis' 'Minor axis' 'Bounding box'};

Get and display spot clusters via kmeans

We initiate to 5 clusters. It could be more or less, and the best number of clusters is to be determined… With clustering we can find spots with a similar size and shape – as well as spots that are around a certain geometrical position, which determines the isoelectric point of the proteins they belong to and their size. Here we care about the size and shape.

[kIdx kCenter] = kmeans(spotProperties, 5);
gelClusters = gelMaskLabel;for i = 1:nSpots
    gelClusters(find(gelMaskLabel == i)) = kIdx(i);
end
figure, imshow(label2rgb(gelClusters, 'jet')), title('Spot Clusters in the gel');

Overlay spot clusters with original image

We use an indexed image of the original gray gel image to overlay with the clusters.

[gelGrayInd gelGrayMap] = gray2ind(gelGray,  128);
gelClustersOverlay = gelGrayInd;

% create 5 more colors in the original map
myExtMap = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 1 0 1; 0 1 1];
theNewMap = [gelGrayMap; myExtMap];

for i = 1:nSpots
    gelClustersOverlay(find(gelMaskLabel == i)) = 128 + kIdx(i);
end
figure, imshow(gelClustersOverlay, theNewMap), title('Spot Clusters overlay');