Commit f6b32a39 authored by Assaf's avatar Assaf

DeBias for CME

parents
function [] = DeBiasCmeCond1Cond2(baseDirname,masterSlaves,params,cond1Str,cond2Str)
close all;
dirnameCond1 = [baseDirname filesep cond1Str];
dirnameCond2 = [baseDirname filesep cond2Str];
if ~exist(dirnameCond1,'dir')
error([dirnameCond1 'does not exists!']);
end
if ~exist(dirnameCond2,'dir')
error([dirnameCond2 'does not exists!']);
end
master = masterSlaves{1};
slave1 = masterSlaves{2};
slave2 = masterSlaves{3};
%% Set parameters
manualAnnotateROIMain(dirnameCond1,params.doAnnotate,master,slave1,slave2);
manualAnnotateROIMain(dirnameCond2,params.doAnnotate,master,slave1,slave2);
close all;
doColocalizationDirectories(dirnameCond1,params,master,slave1,slave2);
doColocalizationDirectories(dirnameCond2,params,master,slave1,slave2);
if params.doMetaAnalysis
metaAnalysisDeBiasCMECond1Cond2(baseDirname,cond1Str,cond2Str,params,master,slave1,slave2);
end
end
%% Meta analysis: 2 conditions
function [] = metaAnalysisDeBiasCMECond1Cond2(baseDirname,cond1Str,cond2Str,params,master,slave1,slave2)
close all;
hasSlave2 = ~isempty(slave2);
addpath(genpath('/home2/azaritsky/code/DeBiasUtils'));
addpath(genpath('/home2/azaritsky/code/extern/export_fig'));
%% Input
dirnameCond1 = [baseDirname filesep cond1Str];
dirnameCond2 = [baseDirname filesep cond2Str];
% MasterSlave1
load([dirnameCond1 filesep 'out.mat']);
GI_Cond1_MasterSlave1 = GIsMasterSlave1(GIsMasterSlave1<7);
LI_Cond1_MasterSlave1 = LIsMasterSlave1(GIsMasterSlave1<7);
CORR_Cond1_MasterSlave1 = CORRsMasterSlave1(GIsMasterSlave1<7);
nCCPs_Cond1_MasterSlave1 = nCCPs;
load([dirnameCond2 filesep 'out.mat']);
GI_Cond2_MasterSlave1 = GIsMasterSlave1(GIsMasterSlave1<7);
LI_Cond2_MasterSlave1 = LIsMasterSlave1(GIsMasterSlave1<7);
CORR_Cond2_MasterSlave1 = CORRsMasterSlave1(GIsMasterSlave1<7);
nCCPs_Cond2_MasterSlave1 = nCCPs;
% MasterSlave2
if hasSlave2
load([dirnameCond1 filesep 'out.mat']);
GI_Cond1_MasterSlave2 = GIsMasterSlave2(GIsMasterSlave2<7);
LI_Cond1_MasterSlave2 = LIsMasterSlave2(GIsMasterSlave2<7);
CORR_Cond1_MasterSlave2 = CORRsMasterSlave2(GIsMasterSlave2<7);
nCCPs_Cond1_MasterSlave2 = nCCPs;
load([dirnameCond2 filesep 'out.mat']);
GI_Cond2_MasterSlave2 = GIsMasterSlave2(GIsMasterSlave2<7);
LI_Cond2_MasterSlave2 = LIsMasterSlave2(GIsMasterSlave2<7);
CORR_Cond2_MasterSlave2 = CORRsMasterSlave2(GIsMasterSlave2<7);
nCCPs_Cond2_MasterSlave2 = nCCPs;
end
%% Output
outdir = baseDirname;
%% LI as function of GI
% MasterSlave1
[rho1_MasterSlave1,pval1_MasterSlave1] = corr(LI_Cond1_MasterSlave1',GI_Cond1_MasterSlave1');
[rho2_MasterSlave1,pval2_MasterSlave1] = corr(LI_Cond2_MasterSlave1',GI_Cond2_MasterSlave1');
fprintf(sprintf('MasterSlave1 LI/GI (n,m,rho,pval): %s(%d,%d,%.2f,%.5f), %s(%d,%d,%.2f,%.5f)\n',...
cond1Str,length(LI_Cond1_MasterSlave1),nCCPs_Cond1_MasterSlave1,rho1_MasterSlave1,pval1_MasterSlave1,...
cond2Str,length(LI_Cond2_MasterSlave1),nCCPs_Cond2_MasterSlave1,rho2_MasterSlave1,pval2_MasterSlave1));
% MasterSlave2
if hasSlave2
[rho1_MasterSlave2,pval1_MasterSlave2] = corr(LI_Cond1_MasterSlave2',GI_Cond1_MasterSlave2');
[rho2_MasterSlave2,pval2_MasterSlave2] = corr(LI_Cond2_MasterSlave2',GI_Cond2_MasterSlave2');
fprintf(sprintf('MasterSlave2 LI/GI (n,m,rho,pval): %s(%d,%d,%.2f,%.5f), %s(%d,%d,%.2f,%.5f)\n',...
cond1Str,length(LI_Cond1_MasterSlave2),nCCPs_Cond1_MasterSlave2,rho1_MasterSlave2,pval1_MasterSlave2,...
cond2Str,length(LI_Cond2_MasterSlave2),nCCPs_Cond2_MasterSlave2,rho2_MasterSlave2,pval2_MasterSlave2));
end
%% Accumulate
labels = [1*ones(length(GI_Cond1_MasterSlave1),1);2*ones(length(GI_Cond2_MasterSlave1),1)];
featsGILI_MasterSlave1 = [[GI_Cond1_MasterSlave1',LI_Cond1_MasterSlave1'];[GI_Cond2_MasterSlave1',LI_Cond2_MasterSlave1']];
[clsGILI_MasterSlave1,resuberrorGILI_MasterSlave1,nErrorsGILI_MasterSlave1,RGILI_MasterSlave1,normRGILI_MasterSlave1] = ...
discreminateColocalization(featsGILI_MasterSlave1,labels,[outdir 'GILI_confusion_MasterSlave1.eps']);
featsCORR_MasterSlave1 = [CORR_Cond1_MasterSlave1';CORR_Cond2_MasterSlave1'];
[clsCORR_MasterSlave1,resuberrorCORR_MasterSlave1,nErrorsCORR_MasterSlave1,RCORR_MasterSlave1,normRCORR_MasterSlave1] = ...
discreminateColocalization(featsCORR_MasterSlave1,labels,[outdir 'CORR_confusion_MasterSlave1.eps']);
if hasSlave2
featsGILI_MasterSlave2 = [[GI_Cond1_MasterSlave2',LI_Cond1_MasterSlave2'];[GI_Cond2_MasterSlave2',LI_Cond2_MasterSlave2']];
[clsGILI_MasterSlave2,resuberrorGILI_MasterSlave2,nErrorsGILI_MasterSlave2,RGILI_MasterSlave2,normRGILI_MasterSlave2] = ...
discreminateColocalization(featsGILI_MasterSlave2,labels,[outdir 'GILI_confusion_MasterSlave2.eps']);
featsCORR_MasterSlave2 = [CORR_Cond1_MasterSlave2';CORR_Cond2_MasterSlave2'];
[clsCORR_MasterSlave2,resuberrorCORR_MasterSlave2,nErrorsCORR_MasterSlave2,RCORR_MasterSlave2,normRCORR_MasterSlave2] = ...
discreminateColocalization(featsCORR_MasterSlave2,labels,[outdir 'CORR_confusion_MasterSlave2.eps']);
end
%% LI & GI
plotGILI(...
GI_Cond1_MasterSlave1,LI_Cond1_MasterSlave1,...
GI_Cond2_MasterSlave1,LI_Cond2_MasterSlave1,...
{cond1Str,cond2Str},...
'LI',... % ylabel
clsGILI_MasterSlave1,...
[outdir 'LIGI_MasterSlave1.eps']...
);
if hasSlave2
plotGILI(...
GI_Cond1_MasterSlave2,LI_Cond1_MasterSlave2,...
GI_Cond2_MasterSlave2,LI_Cond2_MasterSlave2,...
{cond1Str,cond2Str},...
'LI',... % ylabel
clsGILI_MasterSlave2,...
[outdir 'LIGI_MasterSlave2.eps']...
);
end
%%
%% statistics
if params.doStstisticalTest
pGILI_MasterSlave1 = permutationTestCLS(featsGILI_MasterSlave1,labels);
pCorr_MasterSlave1 = permutationTestCLS(featsCORR_MasterSlave1,labels);
fprintf(sprintf('MasterSlave1: pval(GI,LI) = %f\n',pGILI_MasterSlave1));
fprintf(sprintf('MasterSlave1: CORR pval = %f\n',pCorr_MasterSlave1));
if hasSlave2
pGILI_MasterSlave2 = permutationTestCLS(featsGILI_MasterSlave2,labels);
pCorr_MasterSlave2 = permutationTestCLS(featsCORR_MasterSlave2,labels);
fprintf(sprintf('MasterSlave2: pval(GI,LI) = %f\n',pGILI_MasterSlave2));
fprintf(sprintf('MasterSlave2: CORR pval = %f\n',pCorr_MasterSlave2));
end
end
end
%% Utility functions
function [cls,resuberror,nErrors,R,normR] = discreminateColocalization(feats,labels,outfile)
cls = fitcdiscr(feats,labels);
resuberror = resubLoss(cls);
nErrors = resuberror * cls.NumObservations;
R = confusionmat(cls.Y,resubPredict(cls));
normR = R ./ repmat(sum(R,2),1,2);
nLabels = length(unique(labels));
fprintf(sprintf('%s: %.2f\n',outfile,resuberror));
h = figure;
imagesc(normR);
hold on;
haxes = findobj(h,'type','axes');
set(haxes,'XTick',1:nLabels);
set(haxes,'YTick',1:nLabels);
% set(haxes,'XTickLabel',kdStrings);
% set(get(gca,'xlabel'),'Rotation',45);
% set(haxes,'YTickLabel',kdStrings);
set(haxes,'FontSize',24);
caxis([0,1]);
colormap jet;
set(h,'Color','w');
hold off;
export_fig(outfile);
end
%%
function [] = plotGILI(...
GI1,LI1,...
GI2,LI2,...
legendsStr,...
ylabelStr,...
clsGILI,...
outFname...
)
minGI = 2; maxGI = 7;
minLI = -1; maxLI = 4;
assert(...
min(GI1) > minGI && max(GI1) < maxGI && ...
min(LI1) > minLI && max(LI1) < maxLI && ...
min(GI2) > minGI && max(GI2) < maxGI && ...
min(LI2) > minLI && max(LI2) < maxLI...
);
colormap jet;
cmap = colormap(hsv(length(legendsStr)));
fontsize = 20;
h = figure;
hold on;
plot(GI1,LI1,'o','MarkerEdgeColor',cmap(1,:),'LineWidth',2,'MarkerSize',8);
plot(GI2,LI2,'o','MarkerEdgeColor',cmap(2,:),'LineWidth',2,'MarkerSize',8);
xlim([minGI,maxGI]);
if strcmp(ylabelStr,'LI')
ylim([minLI,maxLI]);
end
legend(legendsStr{1},legendsStr{2},'FontSize',24);
xlabel('GI','FontSize',fontsize);
ylabel(ylabelStr,'FontSize',fontsize);
haxes = get(h,'CurrentAxes');
set(h,'Color','none');
set(haxes,'FontSize',fontsize);
export_fig(outFname);
hold off;
end
\ No newline at end of file
function [] = DeBiasCmeTemplate()
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'));
addpath(genpath('/home2/azaritsky/code/common/detectionAlgorithms'));
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
% Update list of experiments
experiments = {...
{'TFRAP2F174','TFR','AP2F174'},... % #1
};
%% Set master / slave channels
prompt = {'Master:','Slave1:','Slave2:'};
dlg_title = 'Set master/slave channels';
num_lines = 1;
defaultans = {'488','561',''};
masterSlaves = inputdlg(prompt,dlg_title,num_lines,defaultans);
%% Set parameters
prompt = {'annotation:','detection:','colocalization:','meta analysis:'};
dlg_title = 'Set parameters';
num_lines = 1;
defaultans = {'0','0','0','1'};
parameters = inputdlg(prompt,dlg_title,num_lines,defaultans);
% params.doAlternativeMeasures = false;
params.doAnnotate = str2double(parameters{1});%true;
params.alwaysDetection = str2double(parameters{2});%false;
params.alwaysColocalization = str2double(parameters{3});%true;
params.intensityFilter = 0;
params.normalizationPrecentile = 5;
params.doMetaAnalysis = str2double(parameters{4});%true;
params.doStstisticalTest = false;
for i = 1 : length(experiments)
tripletDirs = experiments{i};
curDirname = [baseDirname filesep tripletDirs{1} filesep];
cond1 = tripletDirs{2};
cond2 = tripletDirs{3};
DeBiasCmeCond1Cond2(curDirname,masterSlaves,params,cond1,cond2);
end
end
\ No newline at end of file
function [cmePits] = detectCmePits(I488,dirname)
scale=1.8;%getGaussianPSFsigmaFromData(I488,'Display',display);
disp(['Constant scale: ' num2str(scale)]);
% pstruct - detections
% mask - ROI
% imgLM - map of local maximax
[pstruct,mask,imgLM,imgLoG]=pointSourceDetection(I488,scale,'FitMixtures', false,'Alpha',0.05); % low alpha --> less detections
detections.y = pstruct.y;
detections.x = pstruct.x;
cmePits.detections = detections;
h = figure('Visible','Off');
imshow(mat2gray(I488));
hold on
plot(detections.x,detections.y,'ro','MarkerSize',4);
axisHandle= findobj(h,'type','axes');
set(h,'Color','w');
set(axisHandle,'box','off','XMinorTick','off','TickDir','out','YMinorTick','off');
hold off
export_fig([dirname filesep 'detections.eps']);
close all;
end
\ No newline at end of file
function [out] = doColocalization(IMaster,ISlave,cmePits,cmeInds,master,slave,dirname,loggerFname,params)
nIter = 40;
treatStr = [master '_' slave];
xMaster_orig = IMaster(sub2ind(size(IMaster),round(cmePits.detections.y),round(cmePits.detections.x)));
xSlave_orig = ISlave(sub2ind(size(ISlave),round(cmePits.detections.y),round(cmePits.detections.x)));
%% Normalization BEFORE ROI filtering
xMaster = xMaster_orig;
xSlave = xSlave_orig;
%% Hack AZ: looking at CMEs within the ROI + filtering dim structures
% xMaster_roi = xMaster_orig(cmeInds.all);
% xSlave_roi = xSlave_orig(cmeInds.all);
%
% filterTH = prctile(xMaster_roi,params.intensityFilter);
% xMaster = xMaster_roi(xMaster_roi >= filterTH);
% xSlave = xSlave_roi(xMaster_roi >= filterTH);
%%
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));
%% ROI filtering (after normalization)
xMasterNorm = xMasterNorm(cmeInds.all);
xSlaveNorm = xSlaveNorm(cmeInds.all);
%% Number of CCPs analysed
nCCPs = length(xMasterNorm);
[Rho,RhoPval] = corr(xMasterNorm',xSlaveNorm');
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']);
plotColocalizationDistribution(xMasterNorm,xSlaveNorm,-1+0.05:0.05:1-0.05,[dirname treatStr '_Observed.eps']);
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.IMaster = IMaster;
out.ISlave = ISlave;
out.cmePits = cmePits;
out.treatStr = treatStr;
out.xMasterNorm = xMasterNorm;
out.xSlaveNorm = xSlaveNorm;
out.meanAlignment = mean(xMasterNorm - xSlaveNorm);
out.jointDistribution = jointDistribution;
out.nCCPs = nCCPs;
fontsize = 24;
FPosition = [0 0 500 500];
APosition = [0.30 0.30 0.45 0.45];
colorBarYTick = [0,0.01];
h = plotJointDistribution(jointDistribution,colorBarYTick,fontsize);
axisHandle= findobj(h,'type','axes');
hColorbar = colorbar;
set(hColorbar,'YTick',colorBarYTick);
set(hColorbar,'YTickLabel',colorBarYTick.*100);
set(h,'Color','w','Position',FPosition,'PaperPositionMode','auto');
set(axisHandle,'Position',APosition,'box','off','XMinorTick','off','TickDir','out','YMinorTick','off','FontSize',fontsize);
set(get(axisHandle,'XLabel'),'FontSize',fontsize); set(get(axisHandle,'YLabel'),'FontSize',fontsize);
export_fig([dirname treatStr '_jointDistribution.eps']);
% Local/Global index
out.GIs = nan(1,nIter);
out.LIs = nan(1,nIter);
out.PIValidation = nan(1,nIter);
for j = 1 : nIter
[GI, PI, PIValidation] = assessGlobalVsLocalColocalization(xMasterNorm,xSlaveNorm);
out.GIs(j) = GI;
out.LIs(j) = PI;
out.PIValidation(j) = PIValidation;
end
out.Rho = Rho;
out.RhoPval = RhoPval;
out.meanGI = mean(out.GIs);
out.meanLI = mean(out.LIs);
outStr = sprintf('%s (n = %d): local = %.1f, global = %.1f\n',[dirname treatStr],nCCPs,out.meanLI,out.meanGI);
fprintf(outStr);
logger = fopen(loggerFname,'a+');
fprintf(logger,outStr);
fclose(logger);
close all;
end
%% Utilities
function h = plotJointDistribution(jointDist,xaxisVals,fontsize)
h = figure('Visible','Off');
imagescnan(jointDist);
% imagesc(jointDist);
hold on;
colormap('jet');
caxis(xaxisVals); %colorbar;
% haxes = get(h,'CurrentAxes');
haxes = findobj(h,'type','axes');
set(haxes,'XTick',[1,26,50]);
set(haxes,'YTick',[1,26,50]);
set(haxes,'XTickLabel',[0,0.5,1]);
set(haxes,'YTickLabel',[1,0.5,0]);
set(haxes,'FontSize',fontsize);
hold off;
end
function [] = plotDistribution(x,measureStr,bins,outFname)
[nelements, xcenters] = hist(x,bins);
xDistribution = nelements ./ sum(nelements);
h = figure('Visible','Off');
hold on;
bar(bins,xDistribution,'r');
xlabel(measureStr,'FontSize',22);
ylabel('Percent','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.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
function [] = plotColocalizationDistribution(x,y,bins,outFname)
diff = x - y;
[nelements, xcenters] = hist(diff,bins);
colocalizationDistribution = nelements ./ sum(nelements);
h = figure('Visible','Off');
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
% Colocalization for one cells + counting ccps
function [outLIGI] = doColocalizationCell(curdir,loggerFname,params,master,slave)
outFname =[curdir filesep 'out_' master '_' slave '.mat'];
if exist(outFname,'file') && ~params.alwaysColocalization
return;
end
ISlave = double(imread([curdir filesep slave '.tif']));
%% Single cell & detection
detectionsFname = [curdir filesep 'detectionData.mat'];
if exist(detectionsFname,'file') && ~params.alwaysDetection
load(detectionsFname);
else
IMaster = double(imread([curdir filesep master '.tif']));
% Get ROI (if exists)
roiFname = [curdir 'roi.mat'];
if exist(roiFname,'file')
load(roiFname);
else
ROI = true(size(IMaster));
end
stats = regionprops(ROI,'BoundingBox');
foundBB = false;
for i = 1 : length(stats)
tmp = stats(i).BoundingBox;
if tmp(3) > 50 && tmp(4) > 50
if foundBB
error('More than a single ROI');
end
bb = tmp;
foundBB = true;
end
end
x0 = round(bb(1));
x1 = x0 + bb(3)-1;
y0 = round(bb(2));
y1 = y0 + bb(4)-1;
ROI = ROI(y0:y1,x0:x1);
IMaster = IMaster(y0:y1,x0:x1);
ISlave = ISlave(y0:y1,x0:x1);
% IMaster is the CLC images and is used for detecttion
cmePits = detectCmePits(IMaster,curdir);
[cmeInds] = getCmeInds(IMaster,cmePits,ROI);
save(detectionsFname,'IMaster','ROI','cmePits','cmeInds');
end
[out] = doColocalization(IMaster,ISlave,cmePits,cmeInds,master,slave,curdir,loggerFname,params);
outLIGI.GI = out.meanGI;
outLIGI.LI = out.meanLI;
outLIGI.CORR = out.Rho;
outLIGI.nCCPs = out.nCCPs;
save(outFname,'out','cmePits','cmeInds','outLIGI');
close all;
end
\ No newline at end of file
function [] = doColocalizationDirectories(dirname,params,master,slave1,slave2)
i = 1;
if exist([dirname filesep 'out.mat'],'file') && ~params.alwaysColocalization
load([dirname filesep 'out.mat']);
else
loggerFname = [dirname filesep 'log.txt'];
logger = fopen(loggerFname,'w+');
fclose(logger);
GIsMasterSlave1 = [];
LIsMasterSlave1 = [];
CORRsMasterSlave1 = [];
GIsMasterSlave2 = [];
LIsMasterSlave2 = [];
CORRsMasterSlave2 = [];
nCCPs = 0;
countExclude = 0;
while(true)
curDir = [dirname filesep num2str(i) filesep];
if ~exist(curDir,'dir')
countExclude = countExclude + 1;
if countExclude == 100
break;
else
i = i + 1;
continue;
end
end
countExclude = 0;
[outLIGIMasterSlave1] = doColocalizationCell(curDir,loggerFname,params,master,slave1);
[outLIGIMasterSlave2] = doColocalizationCell(curDir,loggerFname,params,master,slave2);
GIsMasterSlave1 = [GIsMasterSlave1 outLIGIMasterSlave1.GI];
LIsMasterSlave1 = [LIsMasterSlave1 outLIGIMasterSlave1.LI];
CORRsMasterSlave1 = [CORRsMasterSlave1 outLIGIMasterSlave1.CORR];
GIsMasterSlave2 = [GIsMasterSlave2 outLIGIMasterSlave2.GI];
LIsMasterSlave2 = [LIsMasterSlave2 outLIGIMasterSlave2.LI];
CORRsMasterSlave2 = [CORRsMasterSlave2 outLIGIMasterSlave2.CORR];
nCCPs = nCCPs + outLIGIMasterSlave1.nCCPs;
i = i + 1;
end
save([dirname filesep 'out.mat'],...
'GIsMasterSlave1','LIsMasterSlave1', 'CORRsMasterSlave1',...
'GIsMasterSlave2','LIsMasterSlave2', 'CORRsMasterSlave2','nCCPs');
end
%% statistical post processing
p_MasterSlave1 = polyfit(GIsMasterSlave1,LIsMasterSlave1,1);
[rho_MasterSlave1, pval_MasterSlave1] = corr(GIsMasterSlave1',LIsMasterSlave1');
[rho_MasterSlave1_corr, pval_MasterSlave1_corr] = corr(GIsMasterSlave1',CORRsMasterSlave1');
%% 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']);
outString = sprintf(...
'Master-Slave1: LI = %.1f (std = %.1f, range: %.1f-%.1f), GI = %.1f (std = %.1f, range: %.1f-%.1f)\n',...
mean(LIsMasterSlave1),std(LIsMasterSlave1),min(LIsMasterSlave1),max(LIsMasterSlave1),...
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));
end
function [] = printPlot(xs,ys,ylabelStr,outFname)
h = figure;
hold on;