Skip to content

Commit 83e5664

Browse files
author
proudot
committed
initial commit
1 parent e72f39c commit 83e5664

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

49 files changed

+4263
-0
lines changed

Diff for: ‎calcDistFromCellEdgeIF.m

+279
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,279 @@
1+
%[out, res] = calcDistFromCellEdgeIF(data, varargin) calculates and plots
2+
% distance and amplitude distributions from the cell edge for spot signals in
3+
% all channels selected in the input
4+
%
5+
% Input:
6+
% data : structure with fields
7+
% .channels : cell array of paths to TIFF files
8+
% .results : path to results file
9+
%
10+
% Options:
11+
% 'Channels' : channel indexes to analyze. Default: [2 3]
12+
% 'BinSize' : histogram bin size, in pixels. Default: 5
13+
%
14+
% Output:
15+
% out : structure with fields
16+
% .dist: distance to cell edge
17+
% .A: amplitude
18+
% .dfeHists: distance from edge histogram
19+
% .ampHists: amplitude histogram (as a function of distance)
20+
% .dfeHistMean_bc: average of bias-corrected dfeHists
21+
% .dfeHistSD_bc: s.d. of bias-corrected dfeHists
22+
% These are the values plotted
23+
% .binc: histogram bin center coordinates
24+
% .ampHistMean_bc: average of bias-corrected ampHists
25+
% .ampHistSD_bc: s.d. of bias-corrected ampHists
26+
%
27+
% Note: 'dist' and 'A' are cell arrays of values for each channel
28+
%
29+
%
30+
% res : structure saved by processFramesIF(), with additional field
31+
% .ex : coordinates (Nx2) of the cell edge
32+
33+
% Francois Aguet, 01/2014
34+
35+
function [out, res] = calcDistFromCellEdgeIF(data, varargin)
36+
37+
ip = inputParser;
38+
ip.CaseSensitive = false;
39+
ip.addRequired('data');
40+
ip.addParameter('BinSize', 5, @isposint);
41+
ip.addParameter('Channels', [2 3]);
42+
ip.addParameter('Display', true, @islogical);
43+
ip.addParameter('Axis', []);
44+
ip.addParameter('Name', '');
45+
ip.addParameter('Hues', []);
46+
ip.addParameter('Names', []);
47+
ip.addParameter('Normalized', true, @islogical);
48+
ip.addParameter('DisplayPlot',true, @islogical);
49+
ip.addParameter('DisplayMode', 'screen');
50+
ip.addParameter('PrintFolder', []);
51+
ip.addParameter('PixelSize', 0.065, @isscalar);
52+
ip.parse(data, varargin{:});
53+
chIdx = ip.Results.Channels;
54+
nc = numel(chIdx);
55+
56+
nd = numel(data);
57+
58+
hues = ip.Results.Hues;
59+
if isempty(hues)
60+
hues = [0.0 0.3];
61+
end
62+
63+
64+
% load data sets
65+
for i = 1:nd
66+
res(i) = load(data(i).results);
67+
end
68+
69+
[ny,nx] = arrayfun(@(i) size(i.mask), res);
70+
xv = 0:ip.Results.BinSize:max([nx ny]/2);
71+
binc = ip.Results.BinSize/2:ip.Results.BinSize:xv(end);
72+
73+
odist = cell(max(chIdx),nd); % distance to cell edge
74+
oA = cell(max(chIdx),nd); % distance to cell edge
75+
odfeHists=cell(max(chIdx),nd);
76+
oampHists=cell(max(chIdx),nd);
77+
78+
parfor i = 1:nd
79+
[oodist,ooA,oodfeHists,ooampHists]=compDist(res(i),i,chIdx,ip.Results,xv);
80+
81+
odist(:,i)=oodist;
82+
oA(:,i)=ooA;
83+
odfeHists(:,i)=oodfeHists;
84+
oampHists(:,i)=ooampHists;
85+
end
86+
out.dist=odist;
87+
out.A=oA;
88+
89+
for c = chIdx
90+
91+
out.dfeHists{c}=vertcat(odfeHists{c,:});
92+
out.ampHists{c}=vertcat(oampHists{c,:});
93+
94+
% correct for observation bias (cell size-dependent), same for all hists/channel
95+
w = sum(out.dfeHists{c}~=0,1);
96+
97+
M = out.dfeHists{c}./repmat(w,[nd 1])*nd;
98+
out.dfeHistsNorm{c}=double(out.dfeHists{c}~=0);
99+
mu = mean(M,1);
100+
s = std(M,[],1)/sqrt(nd);
101+
nidx = isnan(mu);
102+
mu(nidx) = [];
103+
s(nidx) = [];
104+
105+
out.dfeHistMean_bc{c} = mu;
106+
out.s{c} = s;
107+
out.binc{c} = binc(~nidx);
108+
109+
% amplitude
110+
M = out.ampHists{c}./repmat(w,[nd 1])*nd;
111+
out.ampHistsNorm{c}=M;
112+
mu = mean(M,1);
113+
s = std(M,[],1)/sqrt(nd);
114+
nidx = isnan(mu);
115+
mu(nidx) = [];
116+
s(nidx) = [];
117+
out.distBins=xv;
118+
out.ampHistMean_bc{c} = mu;
119+
out.ampHistSD_bc{c} = s;
120+
121+
end
122+
123+
%% Save excel file
124+
125+
126+
if ip.Results.Display
127+
% #objects/cell as a function of distance from cell edge
128+
setupFigure('DisplayMode', ip.Results.DisplayMode, 'Name', ip.Results.Name);
129+
if(~ip.Results.DisplayPlot)
130+
set(gcf(),'Visible','off')
131+
end
132+
133+
hp = zeros(1,nc);
134+
for c = 1:nc
135+
ch = chIdx(c);
136+
xi = out.binc{ch}*ip.Results.PixelSize;
137+
mu = out.dfeHistMean_bc{ch};
138+
s = out.s{ch};
139+
140+
fill([xi xi(end:-1:1)], [mu+s mu(end:-1:1)-s(end:-1:1)], hsv2rgb([hues(c) 0.5 1]),...
141+
'EdgeColor', 'none');
142+
hp(c) = plot(xi, mu, 'Color', hsv2rgb([hues(c) 1 0.9]), 'LineWidth', 1.5);
143+
end
144+
145+
xlabel(['Distance from cell edge (' char(181) 'm)']);
146+
if ip.Results.Normalized
147+
ylabel('Frequency (normalized)');
148+
else
149+
ylabel('# objects/cell');
150+
end
151+
152+
if ~isempty(ip.Results.Axis)
153+
axis(ip.Results.Axis{1});
154+
end
155+
if ~isempty(ip.Results.Names)
156+
hl = legend(hp, ip.Results.Names{:}, 'Location', 'NorthEast');
157+
set(hl, 'Box', 'off');
158+
end
159+
if(~isempty(ip.Results.PrintFolder))
160+
printPNGEPSFIG(gcf(),[ip.Results.PrintFolder filesep], 'AvgObjectCounts')
161+
end
162+
163+
% Average spot amplitudes as a function of distance from cell edge
164+
setupFigure('DisplayMode', ip.Results.DisplayMode, 'Name', ip.Results.Name);
165+
if(~ip.Results.DisplayPlot)
166+
set(gcf(),'Visible','off')
167+
end
168+
for c = 1:nc
169+
ch = chIdx(c);
170+
xi = out.binc{ch}*ip.Results.PixelSize;
171+
mu = out.ampHistMean_bc{ch};
172+
s = out.ampHistSD_bc{ch};
173+
fill([xi xi(end:-1:1)], [mu+s mu(end:-1:1)-s(end:-1:1)], hsv2rgb([hues(c) 0.5 1]),...
174+
'EdgeColor', 'none');
175+
hp(c) = plot(xi, mu, 'Color', hsv2rgb([hues(c) 1 0.9]), 'LineWidth', 1.5);
176+
end
177+
xlabel(['Distance from cell edge (' char(181) 'm)']);
178+
ylabel('Average intensity (A.U.)');
179+
if ~isempty(ip.Results.Axis)
180+
axis(ip.Results.Axis{2});
181+
end
182+
if ~isempty(ip.Results.Names)
183+
hl = legend(hp, ip.Results.Names{:}, 'Location', 'NorthEast');
184+
set(hl, 'Box', 'off');
185+
end
186+
if(~isempty(ip.Results.PrintFolder))
187+
printPNGEPSFIG(gcf(),[ip.Results.PrintFolder filesep], 'AvgIntensities')
188+
end
189+
end
190+
191+
192+
function [oodist,ooA,oodfeHists,ooampHists]=compDist(resData,dataIdx,chIdx,param,xv)
193+
se = strel('disk', 1);
194+
195+
i=dataIdx;
196+
nc=length(chIdx);
197+
198+
[ny, nx] = size(resData.mask);
199+
200+
% cell boundary (including image boundary)
201+
B = bwboundaries(resData.mask);
202+
B = vertcat(B{:}); % [x y] coordinates
203+
bmask = false(ny,nx);
204+
bmask(sub2ind([ny nx], B(:,1), B(:,2))) = true;
205+
206+
% use only detected cell edge (as opposed to image boundary)
207+
% for distance calculation
208+
% emask: cell edge
209+
% bmask: cell edge in image boundary
210+
borderIdx = [1:ny 2*ny:ny:(nx-1)*ny nx*ny:-1:(nx-1)*ny+1 (nx-2)*ny+1:-ny:ny+1];
211+
borderMask = false(ny,nx);
212+
borderMask(borderIdx) = true;
213+
emask = bmask;
214+
emask(borderIdx) = false;
215+
emask = bwmorph(emask, 'skel');
216+
emask = emask | (imdilate(emask,se) & borderMask); % add edge pixel(s) at border
217+
bmask = bmask & ~emask;
218+
219+
[by,bx] = find(bmask);
220+
[ey,ex] = find(emask);
221+
resData.ex = [ex ey];
222+
223+
% distance transform of mask (for normalization)
224+
D = bwdist(~resData.mask);
225+
D2 = bwdist(~padarray(resData.mask, [1 1],0));
226+
D2 = D2(2:end-1,2:end-1);
227+
D(D~=D2) = 0;
228+
D = D-1;
229+
D(D<0) = 0;
230+
231+
% relative probability of observing specific distances from edge
232+
tmp = D(D~=0);
233+
dfeNorm = histc(tmp(tmp<=xv(end)), xv);
234+
dfeNorm = dfeNorm(1:end-1);
235+
236+
if param.Normalized
237+
dfeNorm = dfeNorm/sum(dfeNorm);
238+
else
239+
dfeNorm = dfeNorm/max(dfeNorm);
240+
end
241+
242+
% dirty trick to make it //izable -- PR
243+
nd=1;
244+
oodist = cell(max(chIdx),nd); % distance to cell edge
245+
ooA = cell(max(chIdx),nd); % distance to cell edge
246+
oodfeHists=cell(max(chIdx),1);
247+
ooampHists=cell(max(chIdx),1);
248+
for c = 1:nc
249+
% query cell boundary (free edge + image boundary)
250+
[idx, dist] = KDTreeBallQuery([ex ey; bx by], [resData.ps(chIdx(c)).x' resData.ps(chIdx(c)).y'], max(nx,ny)/2);
251+
252+
% index of points with match;
253+
% results are sorted by increasing distance -> retain first point only
254+
edgeMatchIdx = ~cellfun(@isempty, idx);
255+
idx = cellfun(@(i) i(1), idx(edgeMatchIdx));
256+
dist = cellfun(@(i) i(1), dist(edgeMatchIdx));
257+
258+
% retain only points closer to cell edge than image boundary
259+
dist = dist(idx<=numel(ex));
260+
261+
oodist{chIdx(c),1} = dist; % distance to cell edge
262+
A = resData.ps(chIdx(c)).A(idx<=numel(ex))';
263+
ooA{chIdx(c),1} = A; % amplitude of retained points
264+
265+
[dfeHist, hidx] = histc(dist, xv);
266+
dfeHist = dfeHist(1:end-1)./ dfeNorm;
267+
268+
dfeHist(isnan(dfeHist)) = 0;
269+
if param.Normalized
270+
dfeHist = dfeHist/sum(dfeHist);
271+
end
272+
oodfeHists{chIdx(c)}(1,:) = dfeHist;
273+
274+
% for each bin, calculate average amplitude
275+
ahist = zeros(size(dfeHist));
276+
tmp = arrayfun(@(i) mean(A(hidx==i)), unique(hidx));
277+
ahist(unique(hidx)) = tmp;
278+
ooampHists{chIdx(c)}(1,:) = ahist;
279+
end

Diff for: ‎exportDistFromCellEdgeIFToPrysm.m

+107
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
function [intTable,intTableNorm,countTable,countTableNorm]=exportDistFromCellEdgeIFToPrysm(filename,edgeDistStats,varargin)
2+
% Export stats for prysm
3+
% usage: TODO
4+
% PR 2018
5+
ip = inputParser;
6+
ip.CaseSensitive = false;
7+
ip.addRequired('filename');
8+
ip.addRequired('edgeDistStats');
9+
ip.addParameter('ChannelNames','ch01', @(x) (ischar(x)||iscell(x)));
10+
ip.addParameter('ConditionNames','cond01', @(x) (ischar(x)||iscell(x)));
11+
ip.parse(filename,edgeDistStats,varargin{:});
12+
p=ip.Results;
13+
14+
if(~iscell(edgeDistStats))
15+
edgeDistStats={edgeDistStats};
16+
end
17+
18+
conditionNames=p.ConditionNames;
19+
if(~iscell(conditionNames))
20+
conditionNames={conditionNames};
21+
end
22+
23+
ChannelNames=p.ChannelNames;
24+
if(~iscell(ChannelNames))
25+
ChannelNames={ChannelNames};
26+
end
27+
28+
if(numel(ChannelNames)~=numel(edgeDistStats{1}.dfeHistsNorm))
29+
ChannelNames=arrayfun(@(c) sprintf('Ch%0d',c),1:numel(edgeDistStats{1}.dfeHistsNorm),'unif',0);
30+
end
31+
32+
if(numel(conditionNames)~=numel(edgeDistStats))
33+
conditionNames=arrayfun(@(c) sprintf('Cond%0d',c),1:numel(numel(edgeDistStats)),'unif',0);
34+
end
35+
36+
nCond=length(edgeDistStats);
37+
nCh=numel(ChannelNames);
38+
39+
AllCondStatTable=cell(4,nCond,nCh);
40+
statNames={'count','normCount','int','normInt'};
41+
42+
% condNameRepCell={};
43+
44+
for cIdx=1:length(edgeDistStats)
45+
out=edgeDistStats{cIdx};
46+
xv=out.distBins;
47+
for chIdx=1:nCh
48+
AllCondStatTable(1,cIdx,:)=createStatTables(xv,out.dfeHists,ChannelNames);
49+
AllCondStatTable(2,cIdx,:)=createStatTables(xv,out.dfeHistsNorm,ChannelNames);
50+
AllCondStatTable(3,cIdx,:)=createStatTables(xv,out.ampHists,ChannelNames);
51+
AllCondStatTable(4,cIdx,:)=createStatTables(xv,out.ampHistsNorm,ChannelNames);
52+
end
53+
AllCondStatTable(:,cIdx,:)=cellfun(@(t) prefixVariable(t,conditionNames{cIdx}),AllCondStatTable(:,cIdx,:),'unif',0);
54+
end
55+
56+
% Write in an excel file, a condition and channel per sheet for easier copy
57+
% pasting.
58+
[folder,name,ext]=fileparts(filename);
59+
warning('off','MATLAB:xlswrite:AddSheet')
60+
for sIdx=1:length(statNames)
61+
sfilename=fullfile(folder,[name '-' statNames{sIdx} ext]);
62+
book = matlab.io.spreadsheet.internal.createWorkbook('xls',[], true);
63+
book.save(sfilename)
64+
for condIdx=1:nCond
65+
for chIdx=1:nCh
66+
writetable(AllCondStatTable{sIdx,condIdx,chIdx},sfilename,'Sheet',[conditionNames{condIdx} '-' ChannelNames{chIdx}]);
67+
end
68+
end
69+
end
70+
71+
function T=prefixVariable(T,pref)
72+
T.Properties.VariableNames(2:end)=cellfun(@(n) [pref '_' n], T.Properties.VariableNames(2:end),'unif',0);
73+
74+
75+
function statTableCell=createStatTables(bins,statCell,cellName)
76+
77+
chIdx=1:numel(statCell);
78+
statTableCell=cell(1,numel(chIdx));
79+
80+
81+
for c = chIdx
82+
%% Create 4 tables: Int, raw int, count, rawCount
83+
chName=cellName{c};
84+
85+
cStateTable=array2table(statCell{c}',...
86+
'VariableNames',arrayfun(@(i) sprintf('%s_Cell_%0d',chName,i),1:size(statCell{c},1),'unif',0));
87+
cStateTable.dist=bins(1:height(cStateTable))';
88+
cStateTable=cStateTable(:,[end 1:(end-1)]);
89+
statTableCell{c}=cStateTable;
90+
91+
% if(isempty(statTable))
92+
% statTable=cStateTable;
93+
% else
94+
% statTable=join(statTable,cStateTable,'Keys','dist');
95+
% end
96+
end
97+
98+
99+
function condTable=joinCondTable(condTableCell,condNames)
100+
for c = 1:numel(condTableCell)
101+
condTableCell{c}.Properties.VariableNames(2:end)=cellfun(@(n) [condNames{c} '_' n], condTableCell{c}.Properties.VariableNames(2:end),'unif',0);
102+
end
103+
condTable=condTableCell{1};
104+
for c = 2:numel(condTableCell)
105+
condTable=innerjoin(condTable,condTableCell{c},'Keys','dist');
106+
end
107+

0 commit comments

Comments
 (0)