Skip to content
Snippets Groups Projects
Commit 0b865bfe authored by Sandro Pereira's avatar Sandro Pereira
Browse files

uploading version 2 o Palm analysis

parent 8544c2dc
No related branches found
No related tags found
No related merge requests found
%% Collect Linux\Slurm metadata
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
disp('Job is running on node:')
[~, node] = system('hostname');
disp(node)
disp('Job is run by user:')
[~, user] = system('whoami');
disp(user)
disp('Current slurm jobs of current user:')
[~, sq] = system(['squeue -u ', user]);
disp(sq)
tic
disp(['Start: ' datestr(now, 'yyyymmdd_HHMMSS')])
disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
%%
defineThreshold = false;
if defineThreshold
ThresholdData = load('/mnt/IrisHCS/Data/PatrycjaMulica/VimentinS100B/PKM_20201027_gfap_s100b_20x_TripA53TGCAstro_20201027_150933_out/data.mat');
ThresholdData = ThresholdData.Data;
Q75_S100Vec = ThresholdData.Q75_S100;
Q75_VimentinVec = ThresholdData.Q75_Vimentin;
figure
histogram(Q75_S100Vec, 1000, 'FaceColor', 'r')
hold on
histogram(Q75_VimentinVec, 1000, 'FaceColor', 'g')
legend({'S100', 'Vimentin'})
end
ThresholdS100 = 250;
ThresholdVimentin = 250;
%% Flags
% delete(gcp('nocreate'))
% myCluster = parcluster('local');
% AvailableWorkers = myCluster.NumWorkers;
% if AvailableWorkers >= 28
% pool = parpool(28)
% else
% pool = parpool(1)
% end
addpath(genpath('/work/projects/lcsb_hcs/Library/hcsforge'))
if ~exist('InPath') % if Inpath is provided via command line, use that one
%InPath = '/work/projects/lcsb_hcs/Data/PatrycjaMulica/GFAPSF100beta/PKM_20200723_training_optimization_20200806_111823_in';
%InPath = '/mnt/AtlasHCS/YokogawaCV8000Standalone/CellPathfinder/Patrycja/PKM_20201027_gfap_s100b_20x_TripA53TGCAstro_20201027_150933/AssayPlate_PerkinElmer_CellCarrierUltra96';
InPath = '/mnt/IrisHCS/Data/PatrycjaMulica/SP_markers_Test_astro20XRep1andRep2_20230705_154918_in';
end
MesPath = ls([InPath, '/*.mes']); MesPath = MesPath(1:end-1); % remove line break
MetaData = f_CV8000_getChannelInfo(InPath, MesPath);
if ~exist('OutPath') % if Outpath is provided via command line, use that one
%OutPath = '/work/projects/lcsb_hcs/Data/PatrycjaMulica/GFAPSF100beta/PKM_20200723_training_optimization_20200806_111823_out';
%OutPath = '/mnt/D/tmp20200211';
OutPath = '/mnt/IrisHCS/Data/PatrycjaMulica/SP_markers_Test_astro20XRep1andRep2_20230705_154918_out';
end
%% Prepare folders
mkdir(OutPath)
PreviewPath = [OutPath, filesep, 'Previews'];
mkdir(PreviewPath)
%% Log
f_LogDependenciesLinux(mfilename, OutPath)
%% Load Metadata
ObjectsAll = {};
%MetaData = f_CV8000_getChannelInfo(InPath, MesPath);
InfoTable = MetaData.InfoTable{:};
Wells = unique(InfoTable.Well);
fieldProgress = 0;
for w = 1:numel(Wells)
WellThis = Wells{w};
InfoTableThisWell = InfoTable(strcmp(InfoTable.Well, WellThis),:);
FieldsThisWell = unique(InfoTableThisWell.Field);
for f = 1:numel(FieldsThisWell)
fieldProgress = fieldProgress + 1;
FieldThis = FieldsThisWell{f};
InfoTableThisField = InfoTableThisWell(strcmp(InfoTableThisWell.Field, FieldsThisWell{f}),:);
ChannelsThisField = unique(InfoTableThisField.Channel);
ImPaths = cell(1, numel(ChannelsThisField));
for c = 1:numel(ChannelsThisField)
ChannelThis = ChannelsThisField{c};
InfoTableThisChannel = InfoTableThisField(strcmp(InfoTableThisField.Channel,ChannelThis),:);
InfoTableThisChannel = sortrows(InfoTableThisChannel, 'Plane', 'ascend');
chThisPaths = cell(numel(ChannelsThisField),1);
for p = 1:height(InfoTableThisChannel)
chThisPaths{p} = InfoTableThisChannel{p, 'file'}{:};
%for t = 1:height()
end
ImPaths{c} = chThisPaths;
MesFile = MetaData.MeasurementSettingFileName;
end
FieldMetaData{fieldProgress} = {ImPaths, MesFile, Wells{w}, FieldsThisWell{f}};
end
end
disp('Debug point')
FieldMetaDataTable = cell2table(vertcat(FieldMetaData{:}))
for i = 1:numel(FieldMetaData)
%parfor i = 1:numel(FieldMetaData)
%parfor i = 1:2
try
i
ch1files = sort(FieldMetaData{i}{1}{1});
ch1Collector = cellfun(@(x) imread(x), ch1files, 'UniformOutput', false);
ch1 = cat(3,ch1Collector{:}); % vol(ch1, 0, 2000) Hoechst
ch2files = FieldMetaData{i}{1}{2};
ch2Collector = cellfun(@(x) imread(x), ch2files, 'UniformOutput', false);
ch2 = cat(3,ch2Collector{:}); % vol(ch2, 0, 800) GFAP deep red
ch3files = FieldMetaData{i}{1}{3};
ch3Collector = cellfun(@(x) imread(x), ch3files, 'UniformOutput', false);
ch3 = cat(3,ch3Collector{:}); % vol(ch3, 0, 2000) Nestin 488
ch4files = FieldMetaData{i}{1}{4};
ch4Collector = cellfun(@(x) imread(x), ch4files, 'UniformOutput', false);
ch4 = cat(3,ch4Collector{:}); % vol(ch4, 0, 3000) CellMaskO
MesFile = FieldMetaData{i}{2};
WellThis = FieldMetaData{i}{3};
FieldThis = FieldMetaData{i}{4};
Objects = f_imageAnalysis(ch1, ch2, ch3, ch4, WellThis, FieldThis, MesFile, PreviewPath, ThresholdS100, ThresholdVimentin);
ObjectsAll{i} = Objects;
catch
continue
end
end
Data = vertcat(ObjectsAll{:});
save([OutPath, filesep, 'data.mat'], 'Data');
writetable(Data, [OutPath, filesep, 'data.csv'])
disp('Script completed successfully')
function [FileList] = f_LogDependenciesLinux(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
FileThis = regexp(FileList{i}, '/.*/(.*)', 'tokens'); FileThis = FileThis{:}{:};
SavePathThisFile = [SavePath, filesep, FileThis];
copyfile(FileList{i}, SavePathThisFile);
end
end
function [Objects] = f_imageAnalysis(ch1B, ch2DR, ch3G, ch4O, WellThis, FieldThis, MesFile, PreviewPath, ThresholdS100, ThresholdVimentin)
%Summary is fine do previews and clean comments
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
Summary = table();
%% segment nuclei
chNuc = max(ch1B, [], 3); % imtool(chNuc, [])
chNucMed = medfilt2(chNuc); % imtool(chNucMed, [])
NucMask = chNucMed > 150; % imtool(NucMask, [])
%% segment cells on CellMask
chCellMask = max(ch4O, [], 3); % imtool(chCellMask, [])
chCellMaskMed = medfilt2(chCellMask); % imtool(chCellMaskMed, [])
CellMask = chCellMaskMed > 150; % imtool(CellMask, [])
%% Features
Objects = table();
Objects.Well = {WellThis};
Objects.Field = {FieldThis};
Objects.NucArea = sum(NucMask(:));
Objects.CellArea = sum(CellMask(:));
Objects.meanDeepRedInCellMask = mean(ch2DR(CellMask));
Objects.meanGreenInCellMask = mean(ch3G(CellMask));
%% Previews
imSize = size(chNuc);
[BarMask, BarCenter] = f_barMask(100, 0.32393102760889064, imSize, imSize(1)-50, 50, 20);
%it(BarMask)
BPreview = f_imoverlayIris(imadjust(max(ch1B,[],3), [0 0.03], [0 1]), BarMask, [1 1 1]);%imtool(BPreview)
DRPreview = f_imoverlayIris(imadjust(max(ch2DR,[],3), [0 0.03], [0 1]), BarMask, [1 1 1]);%imtool(DRPreview)
GPreview = f_imoverlayIris(imadjust(max(ch3G,[],3), [0 0.09], [0 1]), BarMask, [1 1 1]);%imtool(GPreview)
OPreview = f_imoverlayIris(imadjust(max(ch4O,[],3), [0 0.03], [0 1]), BarMask, [1 1 1]);%imtool(OPreview)
NucPreview = f_imoverlayIris(imadjust(chNuc, [0 0.01], [0 1]), imdilate(bwperim(NucMask),strel('disk', 1)), [0 0 1]);
NucPreview = f_imoverlayIris(NucPreview, BarMask, [1 1 1]);
%imtool(NucPreview)
CellMaskPreview = f_imoverlayIris(imadjust(chCellMask, [0 0.01], [0 1]), imdilate(bwperim(CellMask),strel('disk', 1)), [1 1 0]);
CellMaskPreview = f_imoverlayIris(CellMaskPreview, BarMask, [1 1 1]);
%imtool(CellMaskPreview)
NucPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_Nuc.png'];
CellMaskPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_CellMask.png'];
BPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_B.png'];
DRPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_DR.png'];
GPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_G.png'];
OPreviewPath = [PreviewPath, filesep, WellThis, '_', FieldThis, '_O.png'];
imwrite(NucPreview, NucPreviewPath)
imwrite(CellMaskPreview, CellMaskPreviewPath)
imwrite(BPreview, BPreviewPath)
imwrite(DRPreview, DRPreviewPath)
imwrite(GPreview, GPreviewPath)
imwrite(OPreview, OPreviewPath)
end
function OverlayRGB = f_imoverlayIris(Im, Mask, RGB)
%IMOVERLAY Adds an overlay mask to a color image without using interactive/graphical Matlab funcctionality.
%
% Recommended for Slurm sbatch runs
% Im has to be two dimensional
% Mask has to be twoo dimensional as well
% RGB is a red green blue vector with a range from 0 to 1 for example [1 0 0]
% Example: Im = im2uint16(rand(100)); Mask = Im > 30000; OverlayRGB = f_imoverlayIris(Im, Mask, [1 0.5 0]); imtool(OverlayRGB)
Mask = logical(Mask);
ImClass = class(Im);
switch ImClass
case 'uint8'
MaxPossible = (2^8)-1;
case 'uint16'
MaxPossible = (2^16)-1;
otherwise
disp('Please input the graytone image as uint8 or as uint16')
return
end
if size(Im, 3) == 1
ImR = Im;
ImG = Im;
ImB = Im;
else
ImR = Im(:,:,1);
ImG = Im(:,:,2);
ImB = Im(:,:,3);
end
ImR(Mask) = RGB(1) * MaxPossible;
ImG(Mask) = RGB(2) * MaxPossible;
ImB(Mask) = RGB(3) * MaxPossible;
OverlayRGB = cat(3, ImR, ImG, ImB);
end
\ No newline at end of file
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment