getFluxZ.m

function Z=getFluxZ(solutionsA, solutionsB)
% getFluxZ
%   Calculates the Z scores between two sets of random flux distributions.
%
%   solutionsA      random solutions for the reference condition (as
%                   generated by randomSampling)
%   solutionsB      random solutions for the test condition (as generated
%                   by randomSampling)
%
%   Z               a vector with Z-scores that tells you for each reaction
%                   how likely it is for its flux to have increased (positive sign)
%                   or decreased (negative sign) in the second condition with
%                   respect to the first.
%
% Usage: Z=getFluxZ(solutionsA, solutionsB)

nRxns=size(solutionsA,1);

%Check that the number of reactions is the same in both cases
if nRxns~=size(solutionsB,1)
    EM='The number of reactions must be the same in solutionsA as in solutionsB';
    dispEM(EM);
end

Z=zeros(nRxns,1);

%Calculate the mean and standard deviation for the two cases
mA=mean(solutionsA,2);
mB=mean(solutionsB,2);

%This can lead to OUT OF MEMORY, so do it in segments of 500 reactions
varA=zeros(size(solutionsA,1),1);
for i=1:500:size(solutionsA,1)
    varA(i:min(i+499,size(solutionsA,1)))=var(solutionsA(i:min(i+499,size(solutionsA,1)),:),0,2);
end
varB=zeros(size(solutionsB,1),1);
for i=1:500:size(solutionsB,1)
    varB(i:min(i+499,size(solutionsB,1)))=var(solutionsB(i:min(i+499,size(solutionsB,1)),:),0,2);
end

%If the mean of both solutions are the same then the Z-score is zero
toCheck=mA~=mB;

%If the variance is zero in both cases, then put a very large or very small
%Z-score for the corresponding reactions
I=find(varA==0 & varB==0 & toCheck==true);
toCheck(I)=false;
J=mA(I)>mB(I);
Z(I(J))=100;
Z(I(~J))=-100;
toCheck=find(toCheck);

for i=1:numel(toCheck)
    Z(toCheck(i))=(mB(toCheck(i))-mA(toCheck(i)))/sqrt(varA(toCheck(i))+varB(toCheck(i)));
end

%Shrink very large values
Z=min(Z,100);
Z=max(Z,-100);
end