README.md

+This code has been tested under Windows 10 version 1903 using Matlab 2019b with all toolboxes installed.
+To run this code example you need to keep the folder structure provided in this repository as depicted below:
+
+CloneClassifier
+|
+|_dataIn
+|
+|_dataOut
+|
+|_src
+
+The image data folder needs to be placed into the folder dataIn
+
+To complete the dependencies, please add the following files to the src folder:
+
+-dirrec.m (can be downloaded from 'https://nl.mathworks.com/matlabcentral/fileexchange/15505-recursive-dir')
+-xmlreadstring.m (can be downloaded from 'https://nl.mathworks.com/matlabcentral/fileexchange/46844-deployanywhere-zip')
+
+To trigger code execution run the main.m script located in src from Matlab
\ No newline at end of file

src/f_CloneClassifier.m

+function [Objects, Preview] = f_CloneClassifier(row, col, ch1Green, ch2Red, ch3Blue, SavePath, ShrinkFactor, Barcode)
+%Classify clones and create previews for a plate overview
+%   Detailed explanation goes here
+
+    %% Image analysis
+
+    % Red pos
+    RedLP = imfilter(ch2Red, fspecial('gaussian', 21, 7), 'symmetric');
+    %it(ch2Red)
+    %it(RedLP)
+    RedPosMask = RedLP > 125;
+    %it(RedPosMask)
+
+    % Green Pos
+    GreenLP = imfilter(ch1Green, fspecial('gaussian', 60, 20), 'symmetric');
+    %it(GreenLP)
+    GreenPosMask = GreenLP > 200;
+    %it(GreenPosMask)
+
+    SamplingZone = RedPosMask | GreenPosMask;
+    %it(SamplingZone)
+    
+    BlueLP = imfilter(ch3Blue, fspecial('gaussian', 60, 20), 'symmetric');
+    BluePosMask = BlueLP > 200;
+
+    %% Features
+
+
+
+    GreenVec = ch2Red(SamplingZone);
+    RedVec = ch1Green(SamplingZone);
+    Objects = table();
+    Objects.Area = sum(SamplingZone(:));
+    Objects.AreaRedPos = sum(RedPosMask(:));
+    Objects.AreaGreenPos = sum(GreenPosMask(:));
+    Objects.AreaBluePos = sum(BluePosMask(:));
+    Objects.GreenMean = mean(GreenVec);
+    Objects.RedMean = mean(RedVec);
+    Objects.Barcode = Barcode;
+    Objects.ROW = row;
+    Objects.COL = col;
+
+
+    if Objects.AreaBluePos / Objects.Area > 0.5
+        Class = 'Blue';
+    else
+        if (Objects.AreaRedPos / Objects.Area > 0.9) & ~(Objects.AreaGreenPos / Objects.Area > 0.1) & (Objects.Area >= 2e4)
+            Class = 'Red';
+        elseif ~(Objects.AreaRedPos / Objects.Area > 0.1) & (Objects.AreaGreenPos / Objects.Area > 0.9) & (Objects.Area >= 2e4)
+            Class = 'Green';
+        elseif (Objects.AreaRedPos / Objects.Area > 0.1) & (Objects.AreaGreenPos / Objects.Area > 0.1) & (Objects.Area >= 2e4)
+            Class = 'RedGreen';
+        else
+            Class = 'Neg';
+        end
+    end
+
+    Objects.Class = {Class};
+
+    %% Previews
+    RGB = cat(3, imadjust(ch2Red, [0 0.01], [0 1]), imadjust(ch1Green, [0 0.01], [0 1]),  imadjust(ch3Blue, [0 0.02], [0 1]));
+    % imtool(RGB)
+    Preview = imoverlay(RGB, imdilate(bwperim(SamplingZone), strel('disk', 10)), [1 1 1]);
+    Preview = insertText(Preview, [1 1], Class, 'FontSize', 200, 'BoxColor', 'black', 'TextColor', 'white');
+    Preview = imresize(Preview, 1/ShrinkFactor);
+    % imtool(Preview)
+
+end
+

src/f_Local_Read_Files_Timeseries_Mosaic_BT2.m

+function [files, info, areaNames, plateLayout, areaMap] = f_Local_Read_Files_Timeseries_Mosaic_BT2(hcsPath, barcode, timepoint)
+%Author: Paul Antony 20141024
+%Loads image data and metadata from a time series barcode from the hcs drive
+%   hcsPath: Full path to the folder containing the meas folder, e.g., 'S:\OperaQEHS\OperaArchiveCol'
+%   barcode: The barcode to be analyzed
+%   timepoint: the counter of the MEAS folder to be analyzed
+%   Example: [files, info, areaNames, plateLayout, areaMap] = f_Local_Read_Files_Timeseries('Y', 'OB 20141022_MDPD1_3Stains_2Controls', 1)
+
+TimepointsStruct = dir([hcsPath, filesep, barcode]);
+TimepointsTable = struct2table(TimepointsStruct);
+TimepointsTableClean = varfun(@(x) regexp(x,'Meas'), TimepointsTable, 'InputVariables', {'name'});
+TimepointsTableFilter = cellfun(@(x) ~isempty(x), TimepointsTableClean.Fun_name, 'UniformOutput', false);
+TimepointsTableFilter = cell2mat(TimepointsTableFilter);
+TimepointsTableClean2 = TimepointsTable(TimepointsTableFilter, TimepointsTable.Properties.VariableNames);
+TimepointFolders = TimepointsTableClean2.name;
+files = dirrec([hcsPath, filesep, barcode, filesep, TimepointFolders{timepoint}], '.flex');
+[info, areaNames, plateLayout, areaMap] = flexinfoMosaic(files);
+
+end
+
+
+
+

src/f_LogDependencies.m

+function [FileList] = f_LogDependencies(Name, SavePath)
+%Collect all upstream scripts and save them in a folder named Dependencies
+%   Name: Function or script name as string
+%   SavePath: 
+%   Example: f_LogDependencies('Workflow_RosellaLC3_BPF_Summary_20170313.m', 'S:\HCS_Platform\Data\JavierJarazo\Rosella\LC3\Analysis_20170323_143443')
+    
+    SavePath = [SavePath filesep 'Dependencies'];
+    mkdir(SavePath);
+    
+    %% Scan recursively %%
+    FileList = matlab.codetools.requiredFilesAndProducts(Name)';
+    %%%%%%%%%%%%%%%%%%%%%%
+    FileCount = size(FileList,1);
+
+    for i = 1:FileCount
+        try % Windows
+            FileThis = regexp(FileList{i}, '\\.*\\(.*)', 'tokens'); FileThis = FileThis{:}{:};
+        catch % Mac
+            FileThis = regexp(FileList{i}, '/.*/(.*)', 'tokens'); FileThis = FileThis{:}{:};
+        end
+        SavePathThisFile = [SavePath, filesep, FileThis];
+        copyfile(FileList{i}, SavePathThisFile);
+    end
+    
+end
+

src/f_load_mosaic_image3D_Sparse.m

+function [XYmosaic, XYmosaicCell, XPositions, YPositions, ResolutionXY] = f_load_mosaic_image3D_Sparse(files, InfoTable, row, column, channel)
+%UNTITLED2 Summary of this function goes here
+%   Detailed explanation goes here
+        if isfield(table2struct(InfoTable), 'row') % For compatibility with data loading before 20161201
+            InfoRowFilter = (InfoTable.row == row) & (InfoTable.column == column);
+        elseif isfield(table2struct(InfoTable), 'Row') % Recommended use
+            InfoRowFilter = (InfoTable.Row == row) & (InfoTable.Column == column);
+        end
+        InfoSelection = InfoTable(InfoRowFilter, InfoTable.Properties.VariableNames);
+        CoordinateStrings = [InfoSelection.PositionX, InfoSelection.PositionY];
+        CoordinatesNumeric = cellfun(@(x) str2double(x), CoordinateStrings);
+        InfoSelection.PositionX = CoordinatesNumeric(:,1);
+        InfoSelection.PositionY = CoordinatesNumeric(:,2);
+
+        InfoSelectionSorted = sortrows(InfoSelection, {'PositionY', 'PositionX'}, 'descend');
+        
+        XPositions = unique(InfoSelectionSorted.PositionX);
+        YPositions = unique(InfoSelectionSorted.PositionY);
+        ResolutionXY = unique(InfoSelectionSorted.ResolutionXY);
+        XYmosaicCell = {};
+        
+        for f = 1:max(InfoSelectionSorted.field(:)) % loop over fields
+            InfoThisField = InfoSelectionSorted((InfoSelectionSorted.field == f), :);
+            cube = readflexcube(InfoThisField.filename{:} , 'PlaneCount', 1);
+            ImageThisField = cube.data(:, :, 1, [channel:InfoThisField.ChannelCount:(InfoThisField.ChannelCount*InfoThisField.PlaneCount)]); %Hoechst
+            % it(ImageThisField);
+            rowThisField = find(YPositions == InfoThisField.PositionY);
+            columnThisField = find(XPositions == InfoThisField.PositionX);
+            XYmosaicCell{rowThisField, columnThisField} = ImageThisField;
+        end
+        
+        EmptyCells = cellfun(@(x) isempty(x), XYmosaicCell);
+        ImSize = size(XYmosaicCell{min(find(~EmptyCells))});
+        XYmosaicCell(EmptyCells) = {zeros(ImSize, 'uint16')};
+        XYmosaic = cell2mat(XYmosaicCell);
+        % it(XYmosaic)
+
+end
+
+

src/flexfileinfoMosaic.m

+function [flexinfo, info] = flexfileinfoMosaic(fname)
+
+flexinfo.filename = fname;
+
+retries = 5;
+while retries > 0, try %#ok<ALIGN>
+    info = imfinfo(fname);
+    retries = 0;
+catch err
+    retries = retries - 1;
+    if retries == 0
+        rethrow(err)
+    end
+    fprintf('Unable to open file %s, %d retries left...\n', fname, retries);
+end, end
+
+xmlstr = info(1).UnknownTags(1).Value;
+xml = xmlreadstring(xmlstr);
+flexinfo.areaname = char(xml.getElementsByTagName('AreaName').item(0).getTextContent());
+flexinfo.barcode = char(xml.getElementsByTagName('Barcode').item(0).getTextContent());
+attr = xml.getElementsByTagName('WellCoordinate').item(0).getAttributes();
+flexinfo.row = str2double(char(attr.getNamedItem('Row').getNodeValue()));
+flexinfo.column = str2double(char(attr.getNamedItem('Col').getNodeValue()));
+flexinfo.meas = str2double(char(xml.getElementsByTagName('OperaMeasNo').item(0).getTextContent()));
+flexinfo.field = str2double(char( xml.getElementsByTagName('Images').item(0).getElementsByTagName('Sublayout').item(0).getTextContent()));
+flexinfo.timestamp = char(xml.getElementsByTagName('DateTime').item(0).getTextContent());
+flexinfo.PositionX = char(xml.getElementsByTagName('PositionX').item(0).getTextContent());
+flexinfo.PositionY = char(xml.getElementsByTagName('PositionY').item(0).getTextContent());
+flexinfo.PositionZ = char(xml.getElementsByTagName('PositionZ').item(0).getTextContent());
+Channels = {};
+for i = 0:xml.getElementsByTagName('Image').getLength-1
+    Channels{i+1} = char(xml.getElementsByTagName('Image').item(i).getElementsByTagName('CameraRef').item(0).getTextContent());
+end
+flexinfo.ChannelCount = length(unique(Channels));
+flexinfo.PlaneCount = xml.getElementsByTagName('Image').getLength / flexinfo.ChannelCount;
+flexinfo.ResolutionXY = eval(char(xml.getElementsByTagName('Image').item(i).getElementsByTagName('ImageResolutionX').item(0).getTextContent()));
+
+
+end
+

src/flexinfoMosaic.m

+function [info, areanames, areas, map] = flexinfoMosaic(files)
+% FLEXINFO extract metadata from flex files
+% [info, areaNames, plateLayout, areaMap] = FLEXFILES(wildcard) 
+% extracts metadata from the flex files matching wildcard
+% flexinfo is a struct with the following fields: 
+% - filename: full path to the flex file
+% - row, column: position of the flex file in the plate
+% - area: numerical identifier for the plate area
+% - areaname: name of the area from the plate layout
+% areaNames is a cell array with areanames indexed by area ids
+% plateLayout is a matrix with the assignment of area ids to well positions
+% areaMap is a containers.Map mapping area names to area ids
+
+%[~, files] = fileattrib(wildcard);
+%info = struct([]);
+%parfor i=1:length(files)
+for i=1:length(files)
+    info(i) = flexfileinfoMosaic(files{i});
+end
+
+[rows, columns] = flexplateinfo(info(1).filename);
+areas = nan(rows, columns);
+
+map = containers.Map();
+
+k = 1;
+
+areanames = cell(1);
+
+for i=1:length(info)
+    row = info(i).row;
+    column = info(i).column;
+    areaname = info(i).areaname;
+    if isKey(map, areaname)
+        area = map(areaname);
+        areas(row, column) = area;
+        info(i).area = area;
+    else
+        areas(row, column) = k;
+        map(areaname) = k;
+        info(i).area = k;
+        areanames{k} = areaname;
+        k = k+1;
+    end
+end
+
+end
+

src/flexplateinfo.m

+function [rows, cols] = flexplateinfo(fname)
+
+info = imfinfo(fname);
+xmlstr = info(1).UnknownTags(1).Value;
+xml = xmlreadstring(xmlstr);
+cols = str2double(char(xml.getElementsByTagName('YSize').item(0).getTextContent()));
+rows = str2double(char(xml.getElementsByTagName('XSize').item(0).getTextContent()));
+
+end
+

src/main.m

+%% Clear Matlab workspace
+
+clear
+clc
+
+%% User inputs per run
+
+InPath = ['..', filesep, 'dataIn'];
+SavePath = ['..', filesep, 'dataOut'];
+Barcode = {'HO_20170921_Subcloning_YEL1'};
+DisplayChannel = 1;
+
+%% Initialize variables and prepare Java environment
+
+mkdir(SavePath)
+ImPlateSummary = {};
+ObjectsAll = {};
+PreviewCell = {};
+
+%% Analysis
+
+for n = 1:size(Barcode, 2)
+
+    Barcode = Barcode{n};
+    Objects = {};
+
+    %% Document script
+
+    SavePathThisBarcode = [SavePath, filesep, Barcode];
+    AnalysisTimeStamp = datestr(now, 'yyyymmdd_HHMMSS');
+    SavePath = [SavePathThisBarcode, filesep, 'Analysis_', AnalysisTimeStamp];
+    mkdir(SavePathThisBarcode)
+    FileNameShort = mfilename;
+    newbackup = sprintf('%s_log.m',[SavePathThisBarcode, filesep, FileNameShort]);
+    FileNameAndLocation = [mfilename('fullpath')];
+    currentfile = strcat(FileNameAndLocation, '.m');
+    copyfile(currentfile,newbackup);
+    Computer = getenv('Computername');
+    Version = ver();
+    save([SavePathThisBarcode, filesep, 'MatlabVersionAndPC.mat'], 'Computer', 'Version')
+    f_LogDependencies(mfilename, SavePathThisBarcode)
+
+    %% Load data
+    [files, info, areaNames, plateLayout, areaMap] = f_Local_Read_Files_Timeseries_Mosaic_BT2(InPath, Barcode, 1);
+
+    %% Organize images
+    % find images for a given mosaic
+
+    InfoTable = struct2table(info);
+    IDs = table([1:size(InfoTable, 1)]');
+    IDs.Properties.VariableNames = {'ID'};
+    InfoTable = [InfoTable, IDs];
+    s = 0
+    PlateImChannels = {};
+
+        for row = unique(InfoTable.row)'
+
+            for column = unique(InfoTable.column)'
+
+                if sum((InfoTable.row == row) .* (InfoTable.column == column)) == 0
+                    continue % Replacing try catch
+                end
+
+                s = s + 1
+                ChannelCount = 3; % more channels needed?
+                XYmosaicCells = cell(1,ChannelCount);
+
+                for channel = 1:ChannelCount
+                   [XYmosaic, XYmosaicCell, XPositions, YPositions, ResolutionXY] = f_load_mosaic_image3D_Sparse(files, InfoTable, row, column, channel);
+                   XYmosaicCells{1,channel} = XYmosaicCell;
+                   XYmosaics{1,channel} = XYmosaic;
+                   PlateImChannels{row, column, channel} = imresize(XYmosaics{1, channel}, 0.075);
+                end
+                
+                ch1Green = cell2mat(XYmosaicCells{1,1});% green
+                ch2Red = cell2mat(XYmosaicCells{1,2});% red
+                ch3Blue = cell2mat(XYmosaicCells{1,3});% blue
+                
+                %% Image analysis goes here
+                [ObjectsThisWell, PreviewThisWell] = f_CloneClassifier(row, column, ch1Green, ch2Red, ch3Blue, SavePathThisBarcode, 1/0.075, Barcode);
+                ObjectsAll{row, column} = ObjectsThisWell;
+                PreviewCell{row, column} = PreviewThisWell;
+                
+            end
+        end
+        
+        
+        %% Collect Data
+        
+        Summary = vertcat(ObjectsAll{:});
+        writetable(Summary, [SavePathThisBarcode, filesep, 'summary.csv'])
+        writetable(Summary, [SavePathThisBarcode, filesep, 'summary.xlsx'])
+        
+        %% Plate previews
+        
+        PlateImChannelsEmpyIdx = cellfun(@(x) isempty(x), PlateImChannels);
+        firstBestImSize = size(PlateImChannels{min(find(~PlateImChannelsEmpyIdx))});
+        PlateImChannels(PlateImChannelsEmpyIdx) = {zeros(firstBestImSize, 'uint16')};
+        
+        NotGoodsizedCells = cellfun(@(x) ~min(size(x) == firstBestImSize), PlateImChannels, 'UniformOutput', false); % Remove images of wrong size (file damage exception)
+        PlateImChannels(cell2mat(NotGoodsizedCells)) = {zeros(firstBestImSize, 'uint16')};
+        
+        plateDisplayImCh1  = cell2mat(PlateImChannels(:,:,1));
+        plateDisplayImCh2  = cell2mat(PlateImChannels(:,:,2));
+        plateDisplayImCh3  = cell2mat(PlateImChannels(:,:,3));
+        
+        plateGridMask = zeros(size(plateDisplayImCh1), 'uint16');
+        plateRows = size(PlateImChannels, 1);
+        plateCols = size(PlateImChannels, 2);
+        
+        plateGridMask([1:firstBestImSize(1):(firstBestImSize(1) * plateRows)],:) = 1;
+        plateGridMask(:, [1:firstBestImSize(2):(firstBestImSize(2) * plateCols)],:) = 1;
+        plateGridMask = imdilate(plateGridMask, strel('disk',3));
+        
+        PlatePreviewCh1 = imoverlay(imadjust(plateDisplayImCh1, [0 0.01], [0 1]), plateGridMask, [0 1 1]); % 520/35
+        imwrite(PlatePreviewCh1, [SavePathThisBarcode, filesep, 'ch1_Green.png'])
+        
+        PlatePreviewCh2 = imoverlay(imadjust(plateDisplayImCh2, [0 0.01], [0 1]), plateGridMask, [0 1 1]); % 600/40
+        imwrite(PlatePreviewCh2, [SavePathThisBarcode, filesep, 'ch2_Red.png'])
+        
+        PlatePreviewCh3 = imoverlay(imadjust(plateDisplayImCh3, [0 0.02], [0 1]), plateGridMask, [0 1 1]); % 450/50
+        imwrite(PlatePreviewCh3, [SavePathThisBarcode, filesep, 'ch3_Blue.png'])
+        
+        %% Classification preview
+        
+        firstBestRGBSize = size(PreviewCell{min(find(~PlateImChannelsEmpyIdx))});
+        PreviewCell(max(PlateImChannelsEmpyIdx,[],3)) = {zeros(firstBestRGBSize, 'uint8')};
+        PlatePreviewClassification = cell2mat(PreviewCell);
+        PlatePreviewClassification = imoverlay(PlatePreviewClassification, plateGridMask, [0 1 1]);
+        imwrite(PlatePreviewClassification, [SavePathThisBarcode, filesep, 'Classes.png'])
+        save([SavePathThisBarcode, filesep, 'Workspace.mat'])
+        
+end
+

src/readflexcube.m

+function [ image, meta ] = readflexcube(fname, varargin)
+% READFLEXCUBE reads XYZC data from a Flex file
+% [ cube, meta ] = READFLEXCUBE(fname)
+% * fname is the filename of the flex file
+% * channels is a numeric array with the indexes of the channels to read
+%   Indices start with 1, if not specified read all channels
+% * meta is the BioFormats metadata store object
+
+parser = inputParser;
+addParamValue(parser, 'Info', []);
+addParamValue(parser, 'PlaneCount', []);
+addParamValue(parser, 'ChannelCount', []);
+addParamValue(parser, 'PixelSizeZ', []);
+parse(parser, varargin{:});
+args = parser.Results;
+
+[~, f] = fileattrib(fname);
+fname = f.Name;
+
+if not(isempty(args.Info))
+    info = args.Info;
+else
+    info = imfinfo(fname);
+end
+
+sizeX = info(1).Width;
+sizeY = info(1).Height;
+sizeZ = 1;
+sizeC = 1;
+
+if not(isempty(args.PlaneCount))
+    sizeZ = args.PlaneCount;
+end
+
+if not(isempty(args.ChannelCount))
+    sizeC = args.ChannelCount;
+end
+
+totalcount = numel(info);
+
+if totalcount ~= sizeZ*sizeC
+    if isempty(args.PlaneCount)
+       sizeZ = floor(totalcount / sizeC);
+    elseif isempty(args.ChannelCount)
+       sizeC = floor(totalcount / sizeZ);
+    end
+end
+
+cube = zeros([sizeY(), sizeX(), sizeZ, sizeC], 'uint16');
+
+getCZTCoords = @(i) [mod(floor(i / sizeC), sizeZ), ...
+                     mod(i, sizeC), ...
+                     floor(floor(i / sizeC) / sizeZ)];
+
+w = warning('off','all');
+tf = Tiff(fname);
+for i=1:totalcount
+   zct = getCZTCoords(i - 1) + 1;
+   cube(:, :, zct(1), zct(2)) = tf.read();
+   if i < totalcount
+       tf.nextDirectory();
+   end
+end
+warning(w);
+
+image.data = cube;
+image.filename = fname;
+
+xmlstr = info(1).UnknownTags(1).Value;
+xml = xmlreadstring(xmlstr);
+
+image.PixelSizeX = str2double(char(xml.getElementsByTagName('ImageResolutionX').item(0).getTextContent()));
+image.PixelSizeY = str2double(char(xml.getElementsByTagName('ImageResolutionY').item(0).getTextContent()));
+
+if not(isempty(args.PixelSizeZ))
+    pixelZ = args.PixelSizeZ;
+else
+    pixelZ = nan;
+end
+
+image.PixelSizeZ = pixelZ;
+
+meta = xml;
+
+end
+