function [aa_P] = substMHPvalues(subst1, subst2);
% substMHPvalues is used to find significant differences between pooled counts of substitutions 
% by the Mantel-Haenszel two-way test.
% Program by Steinar Thorvaldsen, steinart@math.uit.no, Dec. 2006. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last changes 19. Dec 2006.

% Use:     [aa_P]=substMHPvalues(subst1, subst2);

% Input
% subst:    20 x 20 x levels x strata with counts of aa.
%           Counts of amino acids from 2 different groups of organisms and
%           K strata

% Output
% aa_P:     20 x levels. Array with statistical P-values for the observed differences.

[M N D K]=size(subst1)
% The function is made flexible conserning input parameters:
if nargin < 2 
    if M == N %can use reverse substitutions as background to compare with
      subst2=zeros(N,M,D,K);
      for k=1:K
        for d=1:D
          subst2(:,:,d,k)=subst1(:,:,d,k)'; %so that forward subst and reverse subst will have same indecis
        end;
      end;
    end;
end;
[M2 N2 D2 K2]=size(subst2);
if M2 ~= M | N2 ~= N | D2 ~= D | K2 ~= K
    error('Both substitution matrices must be same size!')
end;

MH2x2xK=zeros(2,2,K);
for d=1:D % Each layer is analysed separately
    for aaNo=1:11
        for k=1:K
            %MH2x2xK(1,:,k) = substPairs2aa(subst1(:,:,d,k), aaNo);
            %MH2x2xK(2,:,k) = substPairs2aa(subst2(:,:,d,k), aaNo); %background
            MH2x2xK(1,:,k) = substPairs2bin(subst1(:,:,d,k), aaNo);
            MH2x2xK(2,:,k) = substPairs2bin(subst2(:,:,d,k), aaNo); %background
        end
        %MH2x2xK
        aa_P(aaNo,d) = MantelHaenTest(MH2x2xK); %the Mantel-Haenszel test
    end
end %layer loop
