function [Pmatrix]=substPairsPvalue2FwRev(subst1, subst2, threshold)
% substPairsPvalue is used to compute P-values for substitutions in
% subst1, compared directly with values in subst2 (obtained from a 
% different dataset).
% Two tailed chi2 test is used with a 2x2-array as input.
%(Fisher exect test or Barnard exact test relevant when low samples)
% Program by Steinar Thorvaldsen, april 2004. Last changes 27. Feb 2005.
% Ref.: http://www.math.uit.no/bi/deltaprot/

% The sign of p-values only used to indicate direction of bias, 
% p positive: over-representation, p negative: under-representation.

% The function is made flexible conserning input parameters:
if nargin < 3
    threshold = 1; % all P-values reported, may be confusing to read..
end

[M N D]=size(subst1);
[M2 N2 D2]=size(subst2);
if M2 ~= M | N2 ~= N | D2 ~= D
    error('Both substitution matrices must be same size!')
end;

Pmatrix=zeros(M,N,D);
for d= 1:D
  for i= 1:M
      for j=1:N
          if i~=j %diagonal element with conserved amino acids
            Pmatrix(i,j,d)=chi2Tests([subst1(i,j,d) subst1(j,i,d); subst2(i,j,d) subst2(j,i,d)]);
            % We use + or - sign just to remember the direction of the relative change:
            if ((subst1(i,j,d)/(subst1(i,j,d)+subst1(j,i,d)))-(subst2(i,j,d)/(subst2(i,j,d)+subst2(j,i,d)))) < 0
            %if ((s1(i,j)/(s1(i,j)+s12))-(s2(i,j)/(s2(i,j)+s22))) < 0
              Pmatrix(i,j,d)=-Pmatrix(i,j,d);
            end        
            if abs(Pmatrix(i,j,d)) > threshold
              Pmatrix(i,j,d) = NaN; %just a value to overlook when the user read the matrix
            end
          end
      end
  end
end         