function subst = substPairsCount (S1, S2, level, gap, indep)
% substPairsCount is used to find counts of substitutions, insertions 
% and deletions in different levels between two aligned sequens groups.
% The direction is from S1 to S2. If S1=S2 the counts are within the group.
% Total sums are appended as an extra level (the overall protein).
% Program by Steinar Thorvaldsen, November 2004. Latest update 7. Oct 2009

% Use: subst = substPairsCount (Swarm, Scold, level, 1, 0);

% Input
% S1, S2:  arrays containing several aligned sequences from two groups of organisms
%          Number x sequencelength x 2
% level:   optional vector containing separating layers of the sequences
%          (like 2D structure, surface internal/external, or structural domains). 
%          Other local levels may be derived directly from the sequences (hydrophobic/charged).
%          The levels will separate the output in different context dependent layers.
% gap:     0 if no gaps are to be recorded (default)
%          1 if gaps (insertions and deletions) are to be recorded.
% indep:   statistical independence of data within each group (both within S1 and S2).
%          1 if assumed independent, here we use accumulated counts from most similar pairs
%          0 if dependences (default), here we use mean counts between S1 and S2.

% Output
% subst: 20 x 20 x layers matrix with the mean or total countings between the sequences.
% with gaps included it will be: 21 x 21 x layers

% The function is made flexible conserning input parameters:
if nargin < 5
    indep = 0; %default. Statistical dependence of sequences within same group
end
if nargin < 4
    gap = 0; %default
end
if gap == 1
    A=21; %use gaps as number 21
else
    A=20; %default
end

% If nessesary, convert alignment from a character to an integer representation:
if ischar (S1)
    S1 = aa2int(S1); 
end
if ischar (S2)
    S2 = aa2int(S2); 
end
% To determine the number of sequences M, and length N:
[M1,N] = size (S1);
[M2,N2] = size (S2);
% To have index vectors for the sequences to compare in the default case (indep=0):
m1=1:M1;
m2=1:M2;
if nargin < 3
    level=ones(1,N,'int8'); % no levels
end
L = size (level,2);
if N ~= N2 | N ~= L
   error('All sequences and input vectors must be same length! Make alignment.')
end;

if (M1==M2 && N==N2) & S1==S2 % particular evaluation to count within one group
    D4=M1*(M2-1); % the number of sequences to compare (not to compare to itself)
    within1=true;
    if indep==1
        error('To count within one group, indep=0 must be selected as input parameter')
    end;
    %(the case M1=M2=1 may be handled here....)
else % to count between two groups
    within1=false;
end
D3=max(level); % the number of layers in the context filter
if D3 > 1
   D3=D3+1; % we then need one extra layer to take care of the sum
end

if indep==1 %In this case we will accumulate counts of seq. pairs with maximum similarity:
        %D4=min(M1,M2)^2;% the number of indep sequence pairs to be compared (between two groups)
        idAll=findSeqSimilarity([S1;S2]);
        idS1S2=idAll(1:M1,M1+1:M1+M2); %only look at identities between groups S1 and S2
        id=idS1S2;
        m1S1=m1;
        m2S2=m2;
        origIndexm1=m1;
        origIndexm2=m2;
        m1=[]; %to be used as index for sequenses to compare
        m2=[];
        [idRow,idCol] = size (id); 
        while min(idRow,idCol)>=1 %for d4=1:D4
            maxSim=max(max(id));
            [row,col] = find (id>=maxSim,1,'first');
            m1=[m1 origIndexm1(row)]; %the resulting vector not sorted in ascending order
            m2=[m2 origIndexm2(col)]; %not sorted
            origIndexm1=setdiff(origIndexm1,origIndexm1(row));
            origIndexm2=setdiff(origIndexm2,origIndexm2(col));
            idNewRows=setdiff(m1S1,m1);
            idNewCols=setdiff(m2S2,m2);
            id=idS1S2(idNewRows,idNewCols); %removes already targeted rows and columns
            [idRow,idCol] = size (id);
        end
        subst3 = zeros (A, A, D3); % seq pair substitutions   (A, A, D3,'int8')
        for i=1:length(m1)
            for n=1:N
              aaNo1 = S1(m1(i),n); 
              aaNo2 = S2(m2(i),n);
              if aaNo1 >0 && aaNo1 <A+1 && aaNo2 >0 && aaNo2 <A+1 && aaNo1 ~= aaNo2 %substitution
                    subst3(aaNo1,aaNo2,level(1,n)) = subst3(aaNo1,aaNo2,level(1,n)) + 1;
              end
            end
        end
        %Appending sums as one extra layer in dim 3 and concatenate array into that layer:
        if D3 > 1
           subst3(:,:,D3)=sum(subst3,3);
        end
        subst = subst3;
    
else % indep=0. In this case we will count all pairs and take the mean value in the end
    %D4=M1*M2; % the number of sequence pairs to compare (between two groups)
    subst4 = zeros (A, A, D3, 1);
    d4=0; % to use as index
    for seq1=m1
        for seq2=m2
            if ~(seq1==seq2 & within1) %skip to compare the sequence with itself
                subst3 = zeros (A, A, D3); % seq pair substitutions   (A, A, D3,'int8')
                d4=d4+1;
                for n=1:N
                    aaNo1 = S1(seq1,n); 
                    aaNo2 = S2(seq2,n);
                    if aaNo1 >0 && aaNo1 <A+1 && aaNo2 >0 && aaNo2 <A+1 && aaNo1 ~= aaNo2 %substitution
                        subst3(aaNo1,aaNo2,level(1,n)) = subst3(aaNo1,aaNo2,level(1,n)) + 1;
                    end
                end
                %Appending sums as one extra layer in dim 3 and concatenate array into that layer:
                if D3 > 1
                    subst3(:,:,D3)=sum(subst3,3);
                end
                subst4(:,:,:,d4) = subst3;
            end
        end
    end
    subst = mean(subst4,4); % mean values in dim 4 of all paired counts
    %subst = round(mean(subst4,4)); % rounding off mean values in dim 4
    %subst = round(median(subst4,4)); % median values in dim 4
end
%Appropriate statistical tests:
%Running chi2Test (when indep=0), or Fisher's exact statistical test (indep=1), on each level:
%pmatrix(:,:,d)=chi2Tests(subst(:,:,d)); %mean counts (statistical dependence)
%pmatrix(:,:,d)=FisherExtest(subst(:,:,d)); %total counts (statistical independence)
%Adjusting significance treshold: FDR...

%For plotting results:
%substPairsPlot(subst);
