diff --git a/matlab/+SSF/drawKymosCondition.m b/matlab/+SSF/drawKymosCondition.m index 45ad1431d0830388a0a22dce1cca806aa14349d6..5b330d04ac011fd0bd3677eae289084dfdac2f16 100644 --- a/matlab/+SSF/drawKymosCondition.m +++ b/matlab/+SSF/drawKymosCondition.m @@ -18,15 +18,22 @@ % % cbLabel = 'ERK SSF'; % targetChannels = 1; -function drawKymosCondition(flist,classList,labelNames,outfile,cbLabel,targetChannels,stride) +function drawKymosCondition(flist,classList,labelNames,outfile,cbLabel,targetChannels,stride,framesPerHour) +if ~exist('xticksLabels','var') + xticksLabels.xticks = []; + xticksLabels.xlabels = {}; +end xClasses = unique(classList); nKeep = 0; clipRange = [2.5,97.5]; +% clipRange = [0.5,99.5]; clipLimits = SSF.getClipLimits(flist,targetChannels,clipRange,nKeep); xyim = []; % one per flist +xticks = []; +xtickLabels = {}; for cc = 1:length(xClasses) idxCC = find(strcmp(classList,xClasses(cc))); flistCC = flist(idxCC); @@ -74,29 +81,43 @@ for cc = 1:length(xClasses) plot([x,x],[y,y+szim(ff,1)],'-k') % don't draw to your right... % plot([x+szim(ff,2),x+szim(ff,2)],[y,y+szim(ff,1)],'-k') - text(x,y,ccLabels{ff},'Interpreter','none','FontSize',6,'BackgroundColor',[0.8,0.8,0.8,0.5],... + text(x,y,ccLabels{ff},'Interpreter','none','FontSize',6,'BackgroundColor',[0.9,0.5,0.7,0.5],... 'VerticalAlignment','top'); + + % label time + tHours = [24:24:szim(ff,2)-12]; + for i = 1:length(tHours) + text(x+tHours(i),y+szim(ff,1),num2str(round(tHours(i)./framesPerHour)),'Interpreter','none','FontSize',6,... + 'VerticalAlignment','bottom'); + end + + + + x = x + szim(ff,2); % date_movie = regexp(flistCC(ff).name,'(\d{4}-\d{2}-\d{2}).+_(\d{2})','tokens'); % text(x,y,[date_movie{1}{1} '_' date_movie{1}{2}],'Interpreter','none','FontSize',6,'BackgroundColor',[0.8,0.8,0.8,0.5]); % plot([x+62,x+62],[y+5,y+size(imq,1)-5],'-r') % treatment = regexp(tblCC.Treatment{ff},'(.*)\w*(','tokens'); % text(x+30,y+20,[date_movie{1}{1} '_' date_movie{1}{2} ' :: ' treatment{1}{1}],'Interpreter','none','FontSize',6,'BackgroundColor',[0.8,0.8,0.8,0.7]); + end - + + % set background (index 0, colormap(1)) to white cm = parula(); cm(1,:) = [1,1,1]; colormap(cm); axis on - set(gca,'xtick',[]); + set(gca,'xtick',xticksLabels.xticks); + set(gca,'XTickLabel',xticksLabels.xlabels); set(gca,'ytick',[]); - ylabel('xy \Downarrow x'); + ylabel('xyz \Downarrow x'); xlabel('time'); fx = strfind(flist(ff).name,'.LEVER'); fx = fx(1)-1; - title(xClasses{cc},'Interpreter','none'); - tick_labels = round(linspace(clipLimits(1),clipLimits(2),5),2); +% title(xClasses{cc},'Interpreter','none'); + tick_labels = round(linspace(clipLimits(1),clipLimits(2),5),2,'significant'); cb = colorbar('Ticks',[1,64,128,192,255],'TickLabels',tick_labels); cb.Label.String = cbLabel; if 1 == cc diff --git a/matlab/+SSF/draw_ssf_ncd.m b/matlab/+SSF/draw_ssf_ncd.m index b4d49c1834c0cb95b1e9468064fa690e719ea440..6112032bae5f4ccaa8ed64bfdecb12f6834f8398 100644 --- a/matlab/+SSF/draw_ssf_ncd.m +++ b/matlab/+SSF/draw_ssf_ncd.m @@ -2,23 +2,20 @@ ROOT = '/g/leverjs/Olivier/Agne/march_2023_20x'; flist = dir(fullfile(ROOT,'*.LEVER')); -className = {}; -for ff=1:length(flist) - cn = regexp(flist(ff).name,'(\d+-\d+-\d+)_.*(\d\d)_ori*','tokens'); - cn = [cn{1}{1} '_' cn{1}{2}]; - className{ff} = classMap(cn); -end +classList = get_ssClasses(flist); +[classUnique,ia,ic]= unique(classList); +idxTrue = ic; % targetClass = {'AKT1_E17K','WT','PIK3CA_E545K'} % idx = find(cellfun(@(x) ~isempty(find(strcmp(targetClass,x))),className)); % A=d(idx,idx); % className = className(idx); A = d; -K=3; - +K=6; +[idx,Y] = Cluster.SpectralCluster(A,K); +[csf_mean, csf_std] = CSF.csf_spatial(Y,idxTrue) -[idx,Y] = Cluster.SpectralCluster(A,K); % Y = tsne(Y,'Algorithm','exact','NumDimensions',3); clf;hold off; % if size(Y,2) == 3 @@ -26,21 +23,17 @@ clf;hold off; % else % plot(Y(:,1),Y(:,2),'.') % end -sym = {'ko','cv','m^','r*','gx','<k','b>','gs','bd'}; +sym = {'ko','cv','m^','r*','gx','gs','bd'}; -classList = unique(className); -hx = zeros(size(classList)); +hx = zeros(size(classUnique)); hold on for ff=1:size(Y,1) - idxSym = find(strcmp(classList,className{ff})); - if K==3 + idxSym = find(strcmp(classUnique,classList{ff})); + h1 = plot3(Y(ff,1),Y(ff,2),Y(ff,3),sym{idxSym}); - else - h1 = plot(Y(ff,1),Y(ff,2),sym{idxSym},'UserData',ff); - end if ~hx(idxSym) hx(idxSym) = h1; end @@ -52,5 +45,8 @@ for ff=1:size(Y,1) % end end set(gcf,'color','w') -xlabel('u1');ylabel('u2'),zlabel('u3') -legend(hx(find(hx)),classList(find(hx)),'interpreter','none'); \ No newline at end of file +xlabel('k1');ylabel('k2'),zlabel('k3') +legend(hx(find(hx)),classUnique(find(hx)),'interpreter','none'); +tstr = ['CSF^6_{Erk,Velocity} = ' mat2str([csf_mean, csf_std],2)]; +title(tstr) + diff --git a/matlab/+SSF/getClipLimits.m b/matlab/+SSF/getClipLimits.m index e1303ec949c38bbe821bfc58e7374c31fabbfe1a..423e83296bbdca01b9575d4777b874e278f748a7 100644 --- a/matlab/+SSF/getClipLimits.m +++ b/matlab/+SSF/getClipLimits.m @@ -1,6 +1,6 @@ % map [-1,1] kymographs to color space [0,255] % -function clipLimits = getClipLimits(flist,targetChannelNumbers,clipRange,nKeep) +function clipLimits = getClipLimits(flist,targetChannelNumbers,clipRange,nKeep,tClip) if ~exist('clipRange','var') clipRange= [2.5,97.5]; @@ -8,13 +8,16 @@ end if ~exist('nKeep','var') nKeep = 0; end +if ~exist('tClip','var') + tClip = []; +end p = ljsStartParallel(); kymoPixels = {}; parfor ff=1:length(flist) kp = {}; for c = 1:length(targetChannelNumbers) strDB_ff = fullfile(flist(ff).folder,flist(ff).name); - im_ff = SSF.loadImage(strDB_ff,targetChannelNumbers(c)); + im_ff = SSF.loadImage(strDB_ff,targetChannelNumbers(c), 0, tClip); if nKeep>0 im_ff(im_ff<0)=0; end diff --git a/matlab/+SSF/go_ssf_ncd.m b/matlab/+SSF/go_ssf_ncd.m index b7e82bec83cf6d5216830922423c1ddea3215ab4..e19aaf493800faf5a62c201cc9b2995adc248fe3 100644 --- a/matlab/+SSF/go_ssf_ncd.m +++ b/matlab/+SSF/go_ssf_ncd.m @@ -4,7 +4,7 @@ % AKT_CHANNEL = 2; % H2B_CHANNEL = 3; % targetChannels = [ERK_CHANNEL,AKT_CHANNEL]; -function d=go_ssf_ncd(kROOT,targetChannels,clipRange, nKeep) +function d=go_ssf_ncd(kROOT,targetChannels,clipRange, nKeep, tClip) tStart = tic(); @@ -15,7 +15,15 @@ end if ~exist('clipRange','var') clipRange = [2.5,97.5]; end -flist = dir(fullfile(kROOT,'*.LEVER')); + +if ~exist('tClip','var') + tClip = []; +end +if isstring(kROOT) || ischar(kROOT) + flist = dir(fullfile(kROOT,'*.LEVER')); +else + flist = kROOT; +end d = zeros(length(flist)); kymoPixels = {}; @@ -24,7 +32,7 @@ p = ljsStartParallel(); if iscell(clipRange) clipLimits = clipRange; else - clipLimits = SSF.getClipLimits(flist,targetChannels,clipRange,nKeep); + clipLimits = SSF.getClipLimits(flist,targetChannels,clipRange,nKeep,tClip); end % cmdList = NCD.dParallelCommandList(ones(length(flist)),p.NumWorkers); @@ -44,7 +52,7 @@ parfor i=1:H end strDB_ff = fullfile(flist(ff).folder, flist(ff).name); if ~strcmp(strDB_ff,str_ff_prev) - im_ff = SSF.loadImage(strDB_ff, targetChannels, nKeep); + im_ff = SSF.loadImage(strDB_ff, targetChannels, nKeep, tClip); if iscell(clipLimits) cx = clipLimits{ff}; else @@ -59,7 +67,7 @@ parfor i=1:H dxx(i,j) = d1; else strDB_gg = fullfile(flist(gg).folder,flist(gg).name); - im_gg = SSF.loadImage(strDB_gg,targetChannels, nKeep); + im_gg = SSF.loadImage(strDB_gg,targetChannels, nKeep, tClip); if iscell(clipLimits) cx = clipLimits{gg}; else @@ -74,7 +82,7 @@ parfor i=1:H dxx(i,j) = min(d1,d2); end tElapsed = toc(t0); - fprintf('%s tElapsed = %0.2f\n',mat2str([ff,gg,dxx(i,j)]),tElapsed); +% fprintf('%s tElapsed = %0.2f\n',mat2str([ff,gg,dxx(i,j)]),tElapsed); end end % convert dxx to distance matrix d diff --git a/matlab/+SSF/loadImage.m b/matlab/+SSF/loadImage.m index 54b3f2db8a2deee7fccaed5a752f16801ef64ed6..827e1c330d99bb20e3b86dd49a3bd69cba401a64 100644 --- a/matlab/+SSF/loadImage.m +++ b/matlab/+SSF/loadImage.m @@ -1,5 +1,5 @@ -function im = loadImage(strDB, channelList, nKeep) +function im = loadImage(strDB, channelList, nKeep, tClip) if ~exist('bDouble','var') bDouble = false; @@ -9,6 +9,11 @@ time = 1; % kymo is single time point for c=1:length(channelList) im(:,:,:,c) = leversc.loadImage(strDB,time,channelList(c)); end +% tClip +if exist('tClip','var') && ~isempty(tClip) + tClip(tClip>size(im,3)) = []; + im = im(:,:,tClip,:); +end if ~exist('nKeep','var') return @@ -20,3 +25,4 @@ if nKeep < 0 im(im>0) = 0; im = abs(im); end + diff --git a/matlab/+Segment/frameSegment.m b/matlab/+Segment/frameSegment.m index 58756230a73dbb0f87c25dd2734cd0cc6728d7bf..79c2bdece3aa27415e347a6391b9f0da59d0b04c 100644 --- a/matlab/+Segment/frameSegment.m +++ b/matlab/+Segment/frameSegment.m @@ -143,7 +143,7 @@ for n=1:max(L(:)) newCell.LoG = LoG.sampleCell(newCell,CONSTANTS,imLoG); % check ecc - v = eig(cov(newCell.pts)); + v = eig(cov(double(newCell.pts))); ecc = max(v)/min(v); if ecc>50 continue