By measuring the average average length of the electrophoresis gel image

Question

Reference Information:

My question relates to the extraction function of an electrophoresis gel (see below). In this gel, DNA is loaded from above and allowed to migrate under a voltage gradient. The gel has sieves, so that smaller molecules migrate further than longer molecules, which leads to the separation of DNA in length. The higher the molecule, the longer it is.

Question:

This image has 9 tracks with a separate DNA source. I'm interested in measuring the average location (y-axis value) of each band. I am really new to image processing, but I know MATLAB and I can hardly handle R. I would really appreciate it if someone could show me how to find the average value for each strip.

Answer 1

Here is my attempt. This requires the gels to be good (i.e. Straight paths and the gel should not rotate), but otherwise they will work in general. Please note that there are two parameters depending on the size of the image that will need to be adjusted in order to make this work on images of different sizes.

%# first size-dependent parameter: should be about 1/4th-1/5th %# of the lane width in pixels. minFilterWidth = 10; %# second size-dependent parameter for filtering the %# lane profiles gaussWidth = 5; %# read the image, normalize to 0...1 img = imread('http://img823.imageshack.us/img823/588/gele.png'); img = rgb2gray(img); img = double(img)/255; %# Otsu thresholding to (roughly) find lanes thMsk = img < graythresh(img); 

 %# count the mask-pixels in each columns. Due to %# lane separation, there will be fewer pixels %# between lanes cts = sum(thMsk,1); 

 %# widen the local minima, so that we get a nice %# separation between lanes ctsEroded = imerode(cts,ones(1,minFilterWidth)); %# use imregionalmin to identify the separation %# between lanes. Invert to get a positive mask laneMsk = ~repmat(imregionalmin(ctsEroded),size(img,1),1); 

Image with tracks to be used for analysis

 %# for each lane, create an averaged profile lblMsk = bwlabel(laneMsk); nLanes = max(lblMsk(:)); profiles = zeros(size(img,1),nLanes); midLane = zeros(1,nLanes); for i = 1:nLanes profiles(:,i) = mean(img.*(lblMsk==i),2); midLane(:,i) = mean(find(lblMsk(1,:)==i)); end %# Gauss-filter the profiles (each column is an %# averaged intensity profile G = exp(-(-gaussWidth*5:gaussWidth*5).^2/(2*gaussWidth^2)); G=G./sum(G); profiles = imfilter(profiles,G','replicate'); %' 

 %# find the minima [~,idx] = min(profiles,[],1); %# plot figure,imshow(img,[]) hold on, plot(midLane,idx,'.r') 

Answer 2

Here's my hit on a simple template for an interactive way to do this:

 % Load image img = imread('gel.png'); img = rgb2gray(img); % Identify lanes imshow(img) [x,y] = ginput; % Invert image img = max(img(:)) - img; % Subtract background [xn,yn] = ginput(1); noise = img((yn-2):(yn+2), (xn-2):(xn+2)); noise = mean(noise(:)); img = img - noise; % Calculate means means = (1:size(img,1)) * double(img(:,round(x))) ./ sum(double(img(:,round(x))), 1); % Plot hold on plot(x, means, 'r.') 

Answer 3

The first thing to do is convert the RGB image to grayscale:

 gr = rgb2gray(imread('gelk.png')); 

Then look at the histogram of the image intensity using imhist . Notice anything funny? Use imcontrast(imshow(gr)) to pull out the contrast adjustment tool. I found that eliminating strange material after the main peak of intensity was useful.

The task of image processing can be divided into several stages.

• Split each lane
• Define ("segment") the range in each strip
• Calculate lane location

Step 1 can be performed “manually” if the track bandwidth is guaranteed. Otherwise, the line detection proposed by Hough Transformation is probably the way out. The documentation for Image Processing Toolbox has a really good tutorial on this topic. My code repeats this lesson with the best options for your image. I spent only a few minutes with them, I am sure that you can improve the results by adjusting the parameters further.

Step 2 can be performed in several ways. The easiest way to use it is the Otsu method for creating grayscale images. This method works by defining a threshold that minimizes variance within the -class or, equivalently, maximizes the variance of the interclass. The Otsu method is present in MATLAB as a function of graythresh . If the Otsu method does not work, you can try the multi-level Otsu or several other methods for determining threshold values ​​based on a histogram .

Step 3 can be performed, as you think, by calculating the average y of the segmented pixels of the strip. This is what my code does, although I limited the validation to only the center column of each strip if there was no separation. I am concerned that the result may not be as good as calculating the strip centroid and using its location.

---

Here is my solution:

 function [locations, lanesBW, lanes, cols] = segmentGel(gr) %%# Detect lane boundaries unsharp = fspecial('unsharp'); %# Sharpening filter I = imfilter(gr,unsharp); %# Apply filter bw = edge(I,'canny',[0.01 0.3],0.5); %# Canny edges with parameters [H,T,R] = hough(bw); %# Hough transform of edges P = houghpeaks(H,20,'threshold',ceil(0.5*max(H(:)))); %# Find peaks of Hough transform lines = houghlines(bw,T,R,P,'FillGap',30,'MinLength',20); %# Use peaks to identify lines %%# Plot detected lines above image, for quality control max_len = 0; imshow(I); hold on; for k = 1:length(lines) xy = [lines(k).point1; lines(k).point2]; plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green'); %# Plot beginnings and ends of lines plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow'); plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red'); %# Determine the endpoints of the longest line segment len = norm(lines(k).point1 - lines(k).point2); if ( len > max_len) max_len = len; end end hold off; %%# Use first endpoint of each line to separate lanes cols = zeros(length(lines),1); for k = 1:length(lines) cols(k) = lines(k).point1(1); end cols = sort(cols); %# The lines are in no particular order lanes = cell(length(cols)-1,1); for k = 2:length(cols) lanes{k-1} = im2double( gr(:,cols(k-1):cols(k)) ); %# im2double for compatibility with greythresh end otsu = cellfun(@graythresh,lanes); %# Calculate threshold for each lane lanesBW = cell(size(lanes)); for k = 1:length(lanes) lanesBW{k} = lanes{k} < otsu(k); %# Apply thresholds end %%# Use segmented bands to determine migration distance locations = zeros(size(lanesBW)); for k = 1:length(lanesBW) width = size(lanesBW{k},2); [y,~] = find(lanesBW{k}(:,round(width/2))); %# Only use center of lane locations(k) = mean(y); end 

I suggest that you carefully study not only each output value, but also the results of each step of the function, before using it for real research. To get really good results, you need to read a little about the Hough transforms, Canny edge detection and the Otsu method, and then adjust the parameters. You may also have to change the way lanes are separated; this code assumes that lines will be detected on both sides of the image.

Answer 4

Let me add another implementation similar to @JohnColby's concept, only without manual user interaction:

 %# read image I = rgb2gray(imread('gele.png')); %# middle position of each lane %# (assuming lanes are somewhat evenly spread and of similar width) x = linspace(1,size(I,2),10); x = round( (x(1:end-1)+x(2:end))./2 ); %# compute the mean value across those columns m = mean(I(:,x)); %# find the y-indices of the mean values [~,idx] = min( bsxfun(@minus, double(I(:,x)), m) ); %# show the result figure(1) imshow(I, 'InitialMagnification',100, 'Border','tight') hold on, plot(x, idx, ... 'Color','r', 'LineStyle','none', 'Marker','.', 'MarkerSize',10) 

and applied on a smaller image: