Commit 5b7c5476 authored by Assaf's avatar Assaf

(1) ROI-related bugs solved (2) integrated ZK's version

parent f6b32a39
This diff is collapsed.
This diff is collapsed.
function [] = DeBiasCmeTemplate()
% make sure to addpath to your local DeBiasCME folder!
addpath(genpath('/home2/azaritsky/code/DeBiasCME')); % Update to your code repository
close all;
% addpath(genpath('/home2/azaritsky/code/DeBiasCME'));
addpath(genpath('/home2/azaritsky/code/DeBiasUtils'));
addpath(genpath('/home2/azaritsky/code/extern/export_fig'));
addpath(genpath('/home2/azaritsky/code/common/mathfun'));
......@@ -9,13 +12,19 @@ addpath(genpath('/home2/azaritsky/code/applications/2dActionRecognition/utils'))
addpath(genpath('/home2/azaritsky/code/applications/monolayer/utils'));
%% Here you should hard code the directory and different experiments
baseDirname = '/work/cellbiology/shared/DeBiasAshley/testData'; % Update to data location
% baseDirname = '/work/cellbiology/shared/DeBiasAshley/testData'; % Update to your data folder
baseDirname = '/project/bioinformatics/Danuser_lab/shared/assaf/DeBias/Zuzana/testSW201604';
% Update list of experiments
% experiments = {...
% {'TFRAP2F174','TFR','AP2F174'},... % #1
% };
experiments = {...
{'TFRAP2F174','TFR','AP2F174'},... % #1
{'EGFR_AP2','WT','AP2'},... % #1
};
%% Also include a string for the master and slave channels for printouts (e.g., CLC, EFGR)??
%% Set master / slave channels
prompt = {'Master:','Slave1:','Slave2:'};
dlg_title = 'Set master/slave channels';
......
......@@ -23,32 +23,36 @@ xSlave = xSlave_orig;
lowPctg = params.normalizationPrecentile;
highPctg = 100 - lowPctg;
pLow = prctile(xMaster,lowPctg); pHigh = prctile(xMaster,highPctg); xMaster(xMaster <= pLow) = pLow; xMaster(xMaster >= pHigh) = pHigh;
pLow = prctile(xSlave,lowPctg); pHigh = prctile(xSlave,highPctg); xSlave(xSlave <= pLow) = pLow; xSlave(xSlave >= pHigh) = pHigh;
% normalization 0 - 1
xMasterNorm = (xMaster - min(xMaster)) ./ (max(xMaster)-min(xMaster));
xSlaveNorm = (xSlave - min(xSlave)) ./ (max(xSlave)-min(xSlave));
%% 201604: first select in ROI, then normalize!!
pLow0 = prctile(xMaster(cmeInds),lowPctg); pHigh0 = prctile(xMaster(cmeInds),highPctg);
pLow1 = prctile(xSlave(cmeInds),lowPctg); pHigh1 = prctile(xSlave(cmeInds),highPctg);
%% ROI filtering (after normalization)
xMasterNorm = xMasterNorm(cmeInds.all);
xSlaveNorm = xSlaveNorm(cmeInds.all);
indsInRange = (xMaster >= pLow0) & (xMaster <= pHigh0) & (xSlave >= pLow1) & (xSlave <= pHigh1);
xMasterNorm = (xMaster - pLow0) ./ (pHigh0 - pLow0);
xSlaveNorm = (xSlave - pLow1) ./ (pHigh1 - pLow1);
%% 201604: ROI filtering (after normalization)
xMasterNorm = xMasterNorm(cmeInds & indsInRange);
xSlaveNorm = xSlaveNorm(cmeInds & indsInRange);
%% Number of CCPs analysed
nCCPs = length(xMasterNorm);
[Rho,RhoPval] = corr(xMasterNorm',xSlaveNorm');
out.bins = 0.0125:0.025:1-0.0125;% 201604
plotDistribution(xMasterNorm,'X',0.025:0.025:1-0.025,[dirname treatStr '_X.eps']);
plotDistribution(xSlaveNorm,'Y',0.025:0.025:1-0.025,[dirname treatStr '_Y.eps']);
x0Dist = plotDistribution(xMasterNorm,'X',out.bins,[dirname treatStr '_X.eps']); % 201604: NEW range
x1Dist = plotDistribution(xSlaveNorm,'Y',out.bins,[dirname treatStr '_Y.eps']);
plotColocalizationDistribution(xMasterNorm,xSlaveNorm,-1+0.05:0.05:1-0.05,[dirname treatStr '_Observed.eps']);
plotColocalizationDistribution(xMasterNorm,xSlaveNorm,-1+0.05:0.05:1-0.05,[dirname treatStr '_Observed.eps']); % small bug in the range?
plotColocalizationDistribution(xMasterNorm(randperm(length(xMasterNorm))),xSlaveNorm,-1+0.05:0.05:1-0.05,[dirname treatStr '_Permuted.eps']);
%0:0.02:1,0:0.02:1
binning = 0:0.05:1;
[jointDistribution] = getScatterQuantification(xMasterNorm,xSlaveNorm,binning,binning);
out.jointBins = 0:0.05:1; % 201604
[jointDistribution] = getScatterQuantification(xMasterNorm,xSlaveNorm,out.jointBins,out.jointBins);
out.IMaster = IMaster;
out.ISlave = ISlave;
......@@ -60,6 +64,10 @@ out.meanAlignment = mean(xMasterNorm - xSlaveNorm);
out.jointDistribution = jointDistribution;
out.nCCPs = nCCPs;
% 201604
out.x0Dist = x0Dist;
out.x1Dist = x1Dist;
fontsize = 24;
FPosition = [0 0 500 500];
APosition = [0.30 0.30 0.45 0.45];
......@@ -120,10 +128,11 @@ set(haxes,'FontSize',fontsize);
hold off;
end
function [] = plotDistribution(x,measureStr,bins,outFname)
% 201604: xDistribution is a new output parameter
function [xDistribution] = plotDistribution(x,measureStr,bins,outFname)
[nelements, xcenters] = hist(x,bins);
xDistribution = nelements ./ sum(nelements);
h = figure('Visible','Off');
h = figure;
hold on;
bar(bins,xDistribution,'r');
xlabel(measureStr,'FontSize',22);
......
......@@ -47,7 +47,7 @@ else
% IMaster is the CLC images and is used for detecttion
cmePits = detectCmePits(IMaster,curdir);
[cmeInds] = getCmeInds(IMaster,cmePits,ROI);
[cmeInds] = getCmeInds(cmePits,ROI);
save(detectionsFname,'IMaster','ROI','cmePits','cmeInds');
end
......@@ -58,6 +58,13 @@ outLIGI.LI = out.meanLI;
outLIGI.CORR = out.Rho;
outLIGI.nCCPs = out.nCCPs;
%% 201604
outLIGI.masterDist = out.x0Dist;
outLIGI.slaveDist = out.x1Dist;
outLIGI.jointDistributionMasterSlave = out.jointDistribution;
outLIGI.bins = out.bins;
outLIGI.jointBins = out.jointBins;
save(outFname,'out','cmePits','cmeInds','outLIGI');
close all;
......
function [] = doColocalizationDirectories(dirname,params,master,slave1,slave2)
i = 1;
hasSlave2 = ~isempty(slave2);
if exist([dirname filesep 'out.mat'],'file') && ~params.alwaysColocalization
load([dirname filesep 'out.mat']);
else
......@@ -16,6 +18,18 @@ else
LIsMasterSlave2 = [];
CORRsMasterSlave2 = [];
%% 201604: distributions
masterDist = {};
slave1Dist = {};
slave2Dist = {};
binsDist = nan;
jointDistMasterSlave1 = {};
jointDistMasterSlave2 = {};
binsJointDist = nan;
nCells = 0;
%%
nCCPs = 0;
countExclude = 0;
while(true)
......@@ -33,15 +47,34 @@ else
countExclude = 0;
[outLIGIMasterSlave1] = doColocalizationCell(curDir,loggerFname,params,master,slave1);
[outLIGIMasterSlave2] = doColocalizationCell(curDir,loggerFname,params,master,slave2);
if hasSlave2
[outLIGIMasterSlave2] = doColocalizationCell(curDir,loggerFname,params,master,slave2);
end
GIsMasterSlave1 = [GIsMasterSlave1 outLIGIMasterSlave1.GI];
LIsMasterSlave1 = [LIsMasterSlave1 outLIGIMasterSlave1.LI];
CORRsMasterSlave1 = [CORRsMasterSlave1 outLIGIMasterSlave1.CORR];
GIsMasterSlave2 = [GIsMasterSlave2 outLIGIMasterSlave2.GI];
LIsMasterSlave2 = [LIsMasterSlave2 outLIGIMasterSlave2.LI];
CORRsMasterSlave2 = [CORRsMasterSlave2 outLIGIMasterSlave2.CORR];
if hasSlave2
GIsMasterSlave2 = [GIsMasterSlave2 outLIGIMasterSlave2.GI];
LIsMasterSlave2 = [LIsMasterSlave2 outLIGIMasterSlave2.LI];
CORRsMasterSlave2 = [CORRsMasterSlave2 outLIGIMasterSlave2.CORR];
end
% 201604
masterDist{nCells+1} = outLIGIMasterSlave1.masterDist;
slave1Dist{nCells+1} = outLIGIMasterSlave1.slaveDist;
jointDistMasterSlave1{nCells+1} = outLIGIMasterSlave1.jointDistributionMasterSlave;
if hasSlave2
slave2Dist{nCells+1} = outLIGIMasterSlave2.slaveDist;
jointDistMasterSlave2{nCells+1} = outLIGIMasterSlave2.jointDistributionMasterSlave;
end
binsDist = outLIGIMasterSlave1.bins;
binsJointDist = outLIGIMasterSlave1.jointBins;
nCells = nCells + 1;
%%
nCCPs = nCCPs + outLIGIMasterSlave1.nCCPs;
......@@ -50,7 +83,11 @@ else
save([dirname filesep 'out.mat'],...
'GIsMasterSlave1','LIsMasterSlave1', 'CORRsMasterSlave1',...
'GIsMasterSlave2','LIsMasterSlave2', 'CORRsMasterSlave2','nCCPs');
'GIsMasterSlave2','LIsMasterSlave2', 'CORRsMasterSlave2',...
'masterDist','slave1Dist','slave2Dist',...
'jointDistMasterSlave1','jointDistMasterSlave2',...
'binsDist','binsJointDist',...
'nCells','nCCPs'); % 201604
end
%% statistical post processing
......@@ -58,12 +95,20 @@ p_MasterSlave1 = polyfit(GIsMasterSlave1,LIsMasterSlave1,1);
[rho_MasterSlave1, pval_MasterSlave1] = corr(GIsMasterSlave1',LIsMasterSlave1');
[rho_MasterSlave1_corr, pval_MasterSlave1_corr] = corr(GIsMasterSlave1',CORRsMasterSlave1');
% 201604
if hasSlave2
p_MasterSlave2 = polyfit(GIsMasterSlave2,LIsMasterSlave2,1);
[rho_MasterSlave2, pval_MasterSlave2] = corr(GIsMasterSlave2',LIsMasterSlave2');
[rho_MasterSlave2_corr, pval_MasterSlave2_corr] = corr(GIsMasterSlave2',CORRsMasterSlave2');
end
%% output:
% GI vs. LI
printPlot(GIsMasterSlave1,LIsMasterSlave1,'Local alignment index',[dirname filesep 'GIvsLI_Slave1.eps']);
printPlot(GIsMasterSlave2,LIsMasterSlave2,'Local alignment index',[dirname filesep 'GIvsLI_Slave2.eps']);
if hasSlave2
printPlot(GIsMasterSlave2,LIsMasterSlave2,'Local alignment index',[dirname filesep 'GIvsLI_Slave2.eps']);
end
outString = sprintf(...
'Master-Slave1: LI = %.1f (std = %.1f, range: %.1f-%.1f), GI = %.1f (std = %.1f, range: %.1f-%.1f)\n',...
......@@ -71,11 +116,13 @@ outString = sprintf(...
mean(GIsMasterSlave1),std(GIsMasterSlave1),min(GIsMasterSlave1),max(GIsMasterSlave1));
fprintf(sprintf('%s: \n %s',dirname,outString));
outString = sprintf(...
'Master-Slave2: LI = %.1f (std = %.1f, range: %.1f-%.1f), GI = %.1f (std = %.1f, range: %.1f-%.1f)\n',...
mean(LIsMasterSlave2),std(LIsMasterSlave2),min(LIsMasterSlave2),max(LIsMasterSlave2),...
mean(GIsMasterSlave2),std(GIsMasterSlave2),min(GIsMasterSlave2),max(GIsMasterSlave2));
fprintf(sprintf('%s: \n %s',dirname,outString));
if hasSlave2
outString = sprintf(...
'Master-Slave2: LI = %.1f (std = %.1f, range: %.1f-%.1f), GI = %.1f (std = %.1f, range: %.1f-%.1f)\n',...
mean(LIsMasterSlave2),std(LIsMasterSlave2),min(LIsMasterSlave2),max(LIsMasterSlave2),...
mean(GIsMasterSlave2),std(GIsMasterSlave2),min(GIsMasterSlave2),max(GIsMasterSlave2));
fprintf(sprintf('%s: \n %s',dirname,outString));
end
end
function [] = printPlot(xs,ys,ylabelStr,outFname)
......
function [x0AccDist,x1AccDist,jointAccDist] = getAccDistributions(x0Dists,x1Dists,jointDists)
assert(length(x0Dists) == length(x1Dists))
n = length(x0Dists);
x0AccDist = []; x1AccDist = []; jointAccDist = [];
accN = 0;
for in = 1 : n
x0Dists1 = x0Dists{in};
x1Dists1 = x1Dists{in};
jointDists1 = jointDists{in};
m = length(x0Dists1);
assert(length(x0Dists1) == m);
for im = 1 : m
if isempty(x0AccDist)
x0AccDist = x0Dists1{im};
x1AccDist = x1Dists1{im};
jointAccDist = jointDists1{im};
continue;
end
x0AccDist = x0AccDist + x0Dists1{im};
x1AccDist = x1AccDist + x1Dists1{im};
jointAccDist = jointAccDist + jointDists1{im};
end
accN = accN + m;
end
x0AccDist = x0AccDist ./ accN;
x1AccDist = x1AccDist ./ accN;
jointAccDist = jointAccDist ./ accN;
end
function [cmeInds] = getCmeInds(I488,cmePits,ROI)
function [cmeInds] = getCmeInds(cmePits,ROI)
cmePitsROI = getCmePitsInROI(cmePits,ROI);
n = length(cmePitsROI.detections.x);
cmeCLC = I488(sub2ind(size(I488),round(cmePitsROI.detections.y),round(cmePitsROI.detections.x)));
cmeInds.all = find(cmeCLC);
% cmeInds.inds0_40 = find(cmeCLC < prctile(cmeCLC,40)); cmeInds.inds0_40 = cmeInds.inds0_40(randperm(length(cmeInds.inds0_40),round(0.2*n)));
% cmeInds.inds0_40 = [cmeInds.inds0_40, find(cmeCLC >= prctile(cmeCLC,40))];
%
% cmeInds.inds20_60 = find(cmeCLC > prctile(cmeCLC,20) & cmeCLC < prctile(cmeCLC,60)); cmeInds.inds20_60 = cmeInds.inds20_60(randperm(length(cmeInds.inds20_60),round(0.2*n)));
% cmeInds.inds20_60 = [cmeInds.inds20_60, find(cmeCLC <= prctile(cmeCLC,20) | cmeCLC >= prctile(cmeCLC,60))];
%
% cmeInds.inds40_80 = find(cmeCLC > prctile(cmeCLC,40) & cmeCLC < prctile(cmeCLC,80)); cmeInds.inds40_80 = cmeInds.inds40_80(randperm(length(cmeInds.inds40_80),round(0.2*n)));
% cmeInds.inds40_80 = [cmeInds.inds40_80, find(cmeCLC <= prctile(cmeCLC,40) | cmeCLC >= prctile(cmeCLC,80))];
%
% cmeInds.inds60_100 = find(cmeCLC > prctile(cmeCLC,60) & cmeCLC < prctile(cmeCLC,100)); cmeInds.inds60_100 = cmeInds.inds60_100(randperm(length(cmeInds.inds60_100),round(0.2*n)));
% cmeInds.inds60_100 = [cmeInds.inds60_100, find(cmeCLC <= prctile(cmeCLC,60))];
end
%%
function [cmePitsROI] = getCmePitsInROI(cmePits,ROI)
n = length(cmePits.detections.x);
cmePitsROI.detections.x = [];
cmePitsROI.detections.y = [];
cmeInds = false(1,n);
for i = 1 : n
curY = round(cmePits.detections.y(i));
curX = round(cmePits.detections.x(i));
if(ROI(curY,curX))
cmePitsROI.detections.x = [cmePitsROI.detections.x curX];
cmePitsROI.detections.y = [cmePitsROI.detections.y curY];
cmeInds(i) = true;
end
end
end
\ No newline at end of file
function [cmeInds] = getCmeInds(I488,cmePits,ROI)
cmePitsROI = getCmePitsInROI(cmePits,ROI);
n = length(cmePitsROI.detections.x);
cmeCLC = I488(sub2ind(size(I488),round(cmePitsROI.detections.y),round(cmePitsROI.detections.x)));
cmeInds.all = find(cmeCLC);
% cmeInds.inds0_40 = find(cmeCLC < prctile(cmeCLC,40)); cmeInds.inds0_40 = cmeInds.inds0_40(randperm(length(cmeInds.inds0_40),round(0.2*n)));
% cmeInds.inds0_40 = [cmeInds.inds0_40, find(cmeCLC >= prctile(cmeCLC,40))];
%
% cmeInds.inds20_60 = find(cmeCLC > prctile(cmeCLC,20) & cmeCLC < prctile(cmeCLC,60)); cmeInds.inds20_60 = cmeInds.inds20_60(randperm(length(cmeInds.inds20_60),round(0.2*n)));
% cmeInds.inds20_60 = [cmeInds.inds20_60, find(cmeCLC <= prctile(cmeCLC,20) | cmeCLC >= prctile(cmeCLC,60))];
%
% cmeInds.inds40_80 = find(cmeCLC > prctile(cmeCLC,40) & cmeCLC < prctile(cmeCLC,80)); cmeInds.inds40_80 = cmeInds.inds40_80(randperm(length(cmeInds.inds40_80),round(0.2*n)));
% cmeInds.inds40_80 = [cmeInds.inds40_80, find(cmeCLC <= prctile(cmeCLC,40) | cmeCLC >= prctile(cmeCLC,80))];
%
% cmeInds.inds60_100 = find(cmeCLC > prctile(cmeCLC,60) & cmeCLC < prctile(cmeCLC,100)); cmeInds.inds60_100 = cmeInds.inds60_100(randperm(length(cmeInds.inds60_100),round(0.2*n)));
% cmeInds.inds60_100 = [cmeInds.inds60_100, find(cmeCLC <= prctile(cmeCLC,60))];
end
%%
function [cmePitsROI] = getCmePitsInROI(cmePits,ROI)
n = length(cmePits.detections.x);
cmePitsROI.detections.x = [];
cmePitsROI.detections.y = [];
for i = 1 : n
curY = round(cmePits.detections.y(i));
curX = round(cmePits.detections.x(i));
if(ROI(curY,curX))
cmePitsROI.detections.x = [cmePitsROI.detections.x curX];
cmePitsROI.detections.y = [cmePitsROI.detections.y curY];
end
end
end
\ No newline at end of file
% NOT USED?
function [] = plotDeBiasCMEColocalizationDistribution(x,y,bins,outFname)
diff = x - y;
[nelements, xcenters] = hist(diff,bins);
colocalizationDistribution = nelements ./ sum(nelements);
h = figure;
hold on;
bar(bins,colocalizationDistribution,'r');
xlabel('Colocalization diff','FontSize',22);
ylabel('Percent','FontSize',22);
haxes = get(h,'CurrentAxes');
set(haxes,'XLim',[-1,1]);
set(haxes,'XTick',-1:1:1);
set(haxes,'XTickLabel',-1:1:1);
set(haxes,'YLim',[0,0.2]);
set(haxes,'YTick',0:0.1:0.2);
set(haxes,'YTickLabel',0:0.1:0.2);
set(haxes,'FontSize',22);
set(h,'Color','none');
hold off;
export_fig(outFname);
end
\ No newline at end of file
function [] = plotDeBiasCMEDistribution(distribution,bins,measureStr,outFname)
assert(max(distribution) < 0.071);
h = figure;
hold on;
% plot(bins,distribution,'o-');
bar(bins,distribution,'r');
xlabel(measureStr,'FontSize',22);
ylabel('Freq.','FontSize',22);
haxes = get(h,'CurrentAxes');
set(haxes,'XLim',[0,1]);
set(haxes,'XTick',0:0.5:1);
set(haxes,'XTickLabel',0:0.5:1);
set(haxes,'YLim',[0,0.071]);
set(haxes,'YTick',0:0.03:0.06);
set(haxes,'YTickLabel',0:0.03:0.06);
set(haxes,'FontSize',22);
set(h,'Color','w');
hold off;
export_fig(outFname);
% saveas(h,outFname,'epsc');
end
\ No newline at end of file
%% Utilities
function h = plotDeBiasCMEJointDistribution(jointDist,bins,xaxisVals,outFname)
fontsize = 22;
nBins = length(bins);
h = figure;
imagesc(jointDist);
% imagesc(jointDist);
hold on;
colormap('jet');
caxis(xaxisVals); %colorbar;
% haxes = get(h,'CurrentAxes');
haxes = findobj(h,'type','axes');
set(haxes,'XTick',[1,1+round(nBins/2),nBins-1]);
set(haxes,'YTick',[1,1+round(nBins/2),nBins-1]);
set(haxes,'XTickLabel',[0,0.5,1]);
set(haxes,'YTickLabel',[1,0.5,0]);
set(haxes,'FontSize',fontsize);
% axis equal;
hColorbar = colorbar;
set(hColorbar,'YTick',xaxisVals);
set(hColorbar,'YTickLabel',xaxisVals.*100);
set(h,'Color','w','PaperPositionMode','auto');
set(haxes,'box','off','XMinorTick','off','TickDir','out','YMinorTick','off','FontSize',fontsize);
set(get(haxes,'XLabel'),'FontSize',fontsize); set(get(haxes,'YLabel'),'FontSize',fontsize);
% set(h,'Color','w','Position',FPosition,'PaperPositionMode','auto');
% set(haxes,'Position',APosition,'box','off','XMinorTick','off','TickDir','out','YMinorTick','off','FontSize',fontsize);
% set(get(haxes,'XLabel'),'FontSize',fontsize); set(get(axisHandle,'YLabel'),'FontSize',fontsize);
hold off;
export_fig(outFname);
end
\ No newline at end of file
function [AUC1,AUC2] = rocClsDeBias(cls1,cls2,feats1,feats2,labels,outFname)
[labelsPred1,scores1] = predict(cls1,feats1);
[labelsPred2, scores2] = predict(cls2,feats2);
[X1,Y1,T1,AUC1] = perfcurve(labels,scores1(:,1),1);
[X2,Y2,T2,AUC2] = perfcurve(labels,scores2(:,1),1);
fontsize = 24;
h = figure;
xlabel('False positive rate','FontSize',fontsize);
ylabel('True positive rate','FontSize',fontsize);
hold on;
plot(X1,Y1,'-','Color',[255,165,0]./255,'LineWidth',2);
plot(X2,Y2,'-','Color','k','LineWidth',2);
haxes = get(h,'CurrentAxes');
axis(haxes);
set(h,'Color','none');
set(haxes,'FontSize',fontsize);
export_fig(outFname);
hold off;
close all;
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment