getFluxZ

 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)

0001 function Z=getFluxZ(solutionsA, solutionsB)
0002 % getFluxZ
0003 %   Calculates the Z scores between two sets of random flux distributions.
0004 %
0005 %   solutionsA      random solutions for the reference condition (as
0006 %                   generated by randomSampling)
0007 %   solutionsB      random solutions for the test condition (as generated
0008 %                   by randomSampling)
0009 %
0010 %   Z               a vector with Z-scores that tells you for each reaction
0011 %                   how likely it is for its flux to have increased (positive sign)
0012 %                   or decreased (negative sign) in the second condition with
0013 %                   respect to the first.
0014 %
0015 %   Usage: Z=getFluxZ(solutionsA, solutionsB)
0016 
0017 nRxns=size(solutionsA,1);
0018 
0019 %Check that the number of reactions is the same in both cases
0020 if nRxns~=size(solutionsB,1)
0021     EM='The number of reactions must be the same in solutionsA as in solutionsB';
0022     dispEM(EM);
0023 end
0024 
0025 Z=zeros(nRxns,1);
0026 
0027 %Calculate the mean and standard deviation for the two cases
0028 mA=mean(solutionsA,2);
0029 mB=mean(solutionsB,2);
0030 
0031 %This can lead to OUT OF MEMORY, so do it in segments of 500 reactions
0032 varA=zeros(size(solutionsA,1),1);
0033 for i=1:500:size(solutionsA,1)
0034     varA(i:min(i+499,size(solutionsA,1)))=var(solutionsA(i:min(i+499,size(solutionsA,1)),:),0,2);
0035 end
0036 varB=zeros(size(solutionsB,1),1);
0037 for i=1:500:size(solutionsB,1)
0038     varB(i:min(i+499,size(solutionsB,1)))=var(solutionsB(i:min(i+499,size(solutionsB,1)),:),0,2);
0039 end
0040 
0041 %If the mean of both solutions are the same then the Z-score is zero
0042 toCheck=mA~=mB;
0043 
0044 %If the variance is zero in both cases, then put a very large or very small
0045 %Z-score for the corresponding reactions
0046 I=find(varA==0 & varB==0 & toCheck==true);
0047 toCheck(I)=false;
0048 J=mA(I)>mB(I);
0049 Z(I(J))=100;
0050 Z(I(~J))=-100;
0051 toCheck=find(toCheck);
0052 
0053 for i=1:numel(toCheck)
0054     Z(toCheck(i))=(mB(toCheck(i))-mA(toCheck(i)))/sqrt(varA(toCheck(i))+varB(toCheck(i)));
0055 end
0056 
0057 %Shrink very large values
0058 Z=min(Z,100);
0059 Z=max(Z,-100);
0060 end