AnalysisRNAPanomicsTissue.m

%% AnalysisRNAPanomicsTissue
clc
clear variables
close all

%% User Settings
pathData = '/Users/tushevg/Desktop/180817-FC-1R-first slicing-2aq';
channelDapi = 1;
channelMap2 = 2;
channelRna = 3;
localMinRadius = 3;
punctaMinIntensity = 100;



%% Loop over each file in the directory
fileList = getFilesFromDir(pathData, 'tif');
fileResult = regexp(pathData, filesep, 'split');
fileResult = regexprep(fileResult{end},'\ ','_');
fileResult = [fileResult,datestr(now,'yyyymmdd_HHMMSS'),'.txt'];
fileResult = fullfile(pathData, fileResult);
fh = fopen(fileResult, 'w');
fprintf(fh, '# path = %s\n', pathData);
fprintf(fh, '# images = %d\n', length(fileList));
fprintf(fh, '# channel DAPI = %d\n', channelDapi);
fprintf(fh, '# channel MAP2 = %d\n', channelMap2);
fprintf(fh, '# channel RNA = %d\n', channelRna);
fprintf(fh, '# local minima radius = %.2f\n', localMinRadius);
fprintf(fh, '# puncta minimum intensity = %.2f\n', punctaMinIntensity);
fprintf(fh, '# image\tpuncta.soma\tpuncta.neuropil\tarea.soma [pixels]\tarea.neuropil [pixels]\n');

f = 1;
for f = 1 : 3 %length(fileList)
    
    
    % read image and split channels
    img = imread(fileList{f});
    
    imgCell = imadd(img(:, :, channelDapi),...
                      img(:, :, channelMap2));
                  
    imgRna = img(:, :, channelRna);
    
    % check if region exists
    [pathName, fileName] = fileparts(fileList{f});
    fileBorder = fullfile(pathName, [fileName,'.txt']);
    
    if exist(fileBorder, 'file') == 2
        
        fr = fopen(fileBorder, 'r');
        fgetl(fr);
        fgetl(fr);
        txt = textscan(fr, '%n %n','delimiter','\t');
        fclose(fr);
        border = [txt{1},txt{2}];
    else
        hfig = figure('color','w');
        himg = imshow(imgCell,[]);
        hline = impoly(himg.Parent, 'Closed', false);
        border = hline.getPosition;
        close(hfig);
        
        
        fw = fopen(fileBorder, 'w');
        fprintf(fw,'# Border\n');
        fprintf(fw,'# y\tx\n');
        fprintf(fw,'%.2f\t%.2f\n',border');
        fclose(fw);
        
    end
                           
    % create mask
    mask = regionToMask(border,...
                        imgCell,...
                        true);
                    
    % find puncta
    [bw_smt, punctaList_smt] = getPuncta(imgRna, mask, localMinRadius, punctaMinIntensity);
    [bw_npl, punctaList_npl] = getPuncta(imgRna, ~mask, localMinRadius, punctaMinIntensity);
    
    fprintf(fh, '%s\t%d\t%d\t%d\t%d\n',...
            fileName,...
            length(punctaList_smt),...
            length(punctaList_npl),...
            sum(mask(:)),...
            sum(~mask(:)));
    
end

fclose(fh);




%% FUNCTIONS
function fileList = getFilesFromDir(pathName, fileExtension)
    filesQuery = [pathName, filesep, '*.', fileExtension];
    dirInfo = dir(filesQuery);
    fileList = cellfun(@(x) {fullfile(pathName, x)}, {dirInfo.name}');
end

function [xEdges, yEdges] = getEdgePoints(x, y, w, h)
    

    closestPoint_start = getClosestPointToRectangle(x(1), y(1), w, h);
    closestPoint_end = getClosestPointToRectangle(x(end), y(end), w, h);

    xEdges = [closestPoint_start(1); x; closestPoint_end(1)];
    yEdges = [closestPoint_start(2); y; closestPoint_end(2)];
    
end


function closestPoint = getClosestPointToRectangle(x, y, w, h)

    edgePoints = [0, y;...
                  w, y;...
                  x, 0;...
                  x, h];
              
    [~, i] = min([x;...
                  (w - x);...
                  y;...
                  (h - y)]);
              
    closestPoint = edgePoints(i, :);          

end

function [xPolygon, yPolygon] = getPolygonPoints(x, y, w, h)

    cornerPoints = [0,0;...
                    0,h;...
                    w,0;...
                    w,h];

    value = distEuclidean(mean(x), mean(y),...
                              cornerPoints(:,1),...
                              cornerPoints(:,2));
                          
    [~,i] = sort(value);
    cornerPoints = cornerPoints(i,:);
    cornerPoints = cornerPoints(1:2,:);
    
    value = distEuclidean(x(1), y(1),...
                          cornerPoints(:,1),...
                          cornerPoints(:,2));
    [~,i] = sort(value);
    cornerPoints = cornerPoints(i,:);
    
    
    xPolygon = [cornerPoints(1, 1);...
                x;...
                cornerPoints(2, 1)];
            
    yPolygon = [cornerPoints(1,2);...
                y;...
                cornerPoints(2, 2)];
   
end

function value = distEuclidean(x,y,m,n)

    value = sqrt((x - m).^2 + (y - n).^2);
    
end

function mask = regionToMask(border, img, showPlots)

    [h, w] = size(img);

    % locate image edge points
    [xBorder, yBorder] = getEdgePoints(border(:,2),...
                                       border(:,1),...
                                       w,...
                                       h);              
                            
    % generate interpolation base
    t = [0; cumsum(diff(xBorder).^2 + diff(yBorder).^2)];
    ti = linspace(0, t(end), 100)';
    
    xArray = pchip(t, xBorder, ti);
    yArray = pchip(t, yBorder, ti);
    
    % create closed polygon
    [xPoly, yPoly] = getPolygonPoints(xArray,...
                                      yArray,...
                                      w,...
                                      h);
    
    % generate mask                              
    mask = poly2mask(yPoly, xPoly, w, h);
    
    % show data
    if showPlots
        
        colorMask = cat(3, ones(size(img),class(img)),...
                           zeros(size(img),class(img)),...
                           zeros(size(img),class(img)));
        
        figure('color','w');
        imshow(img,[]);
        hold on;
        hmask = imshow(colorMask,[]);
        plot(yArray, xArray,'-','color','g');
        plot(yBorder, xBorder, 'ro');
        hold off;
        set(hmask, 'AlphaData', mask .*0.5);
        
    end

end

function [bw, punctaList] = getPuncta(img, mask, radius, inty)
    
    img = imgaussfilt(img, 0.2);
    bw = (img == imdilate(img, strel('disk', radius))) & ...
                (img >= inty) & ...
                mask;
    rp = regionprops(bw, img, 'WeightedCentroid');
    punctaList = cat(1, rp.WeightedCentroid);
    
end