function [Pmatrix]=substPairsPvalues(subst1, subst2, threshold)
% substPairsPvalue is used to compute P-values for substitutions in
% subst1, with values in subst2 as backgroun (reference).
% 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 31. Oct 2006.
% Ref.: http://www.math.uit.no/bi/deltaprot/

% This procedure has its implicit direction of change:
% FROM index 1 (row) TO index 2 (column).
% The sign of p-values only used to indicate the direction of the 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);
if nargin < 2 
    if M == N %can use reverse substitutions as background to compare with
      subst2=zeros(N,M,D);
      for d=1:D
        subst2(:,:,d)=subst1(:,:,d)'; %so that forward subst and reverse subst will have same indecis
      end
    else %pooled data vector
      error('Pooled substitution data: must have two input vectors!')
    end;
else
    [M2 N2 D2]=size(subst2);
    if M2 == N & N2 == M & D2 == D %a pooled data case with forward and reverse substitutions
      temp2=zeros(M,N,D);
      for d=1:D
        temp2(:,:,d)=subst2(:,:,d)'; 
      end
      subst2=temp2; %so that forward and reverse vectors will have the same indecis
    elseif M2 ~= M | N2 ~= N | D2 ~= D
        error('Both substitution matrices must be same size!')
    end;
end

Pmatrix=zeros(M,N,D);
for d= 1:D
  s1=subst1(:,:,d);
  s2=subst2(:,:,d);
  if M == N % To be sure that all diagonal counts (non-substitutions) are removed:
      s1=s1-diag(diag(s1));
      s2=s2-diag(diag(s2));
  end
  marg_row1=sum(s1,2);
  marg_row2=sum(s2,2);
  
  if N == 1 % Input are column-vectors (with pooled data)
     % Replicates to just to avoid index overflow problems:
     temp1=sum(marg_row1,1);
     temp2=sum(marg_row2,1);
     marg_row1(1:M,N)=temp1; 
     marg_row2(1:M,N)=temp2;
  end
 
  for i= 1:M
      for j=1:N
          if i==j %diagonal element with conserved amino acids
            s12=marg_row1(i);
            s22=marg_row2(i);
            %%Pmatrix(i,j,d)=NaN; %FisherExtest([s1(i,i) marg_row1(i);s2(i,i) marg_row2(i)]);
          else %substututions
            s12=marg_row1(i)-s1(i,j);
            s22=marg_row2(i)-s2(i,j);
            %Pmatrix(i,j,d)=FisherExtest([s1(i,j) marg_row1(i)-s1(i,j); s2(i,j) marg_row2(i)-s2(i,j)]);
          end
          % avoids negative values in the chi2Test caused by round off errors:
          if abs(s12)<eps(single(1)) s12=0; end;
          if abs(s22)<eps(single(1)) s22=0; end;
          Pmatrix(i,j,d)=chi2Tests([s1(i,j) s12; s2(i,j) s22]);
          % We use + or - sign just to remember the direction of the relative change:
          denum1=s1(i,j)+s12;
          denum2=s2(i,j)+s22;
          if denum1 ~= 0 & denum2 ~= 0
             if ((s1(i,j)/denum1)-(s2(i,j)/denum2)) < 0
                Pmatrix(i,j,d)=-Pmatrix(i,j,d);
             end
          end
          if abs(Pmatrix(i,j,d)) > threshold
            %Pmatrix(i,j,d) = NaN; %just a value to easily overlook when the user read the matrix
          end
      end
  end
end
