function aaDfreq = aaDeltaFreq(Sgrouped, level);
% aaDeltaFreq is used to compute the change in mean aa-frequences between 
% proteins in the populations, where the firste group is treated as baseline.
% Program by Steinar Thorvaldsen, april 2004. Last changes 24. Nov 2005.

% Use:      
%           aaDfreq = aaDeltaFreq(Sgrouped, level)
% Input
% Sgrouped: cell array containing the gene sequence from 2 or more different groups of organisms
% level:    a vector containin information on 2D or 3D about the aligned sequences. 
%           The number of layers may be 1 as for unaligned sequences.
% Output
% aaDfreq:  (Ngroups-1) x 20 x levels, a level-separated matrix with aa countings from the N groups.
%           Gaps are neglected.

% Reconstructing the dimension of the data:
Ngroups = size(Sgrouped,1); % the number of groupes
for n=1:Ngroups
    nGroup(n)=size(Sgrouped{n},1); % the number of sequences in each group
end
M = sum(nGroup); % the total number of sequences
[M1,N1] = size(cell2mat(Sgrouped(1))); %N1 must be same in all seq groups

% The function is made flexible conserning input parameters:
if nargin < 2
    level=ones(N1,1,'int8'); % no levels
    D3=1;
% Some alignments may have few (~3) levels (no gaps). May cause problems***:
else % the number of levels should stay constant throught all seq analysis
    global maxlevel;
    D3=maxlevel;
end

% Starting to compute frequences of amino acids:
aa_total = aaCount(cell2mat(Sgrouped), level); % aa compositions for each seq from function aa_count
if D3 > 1
    aa_total=cat(3,aa_total,sum(aa_total,3)); %Concatenate two arrays by appending sums as extra layer
    D3=D3+1;
end
% The frequency change for N groups relative to the baseline group separaten in D3 layers:
aaDfreq = zeros (Ngroups-1,20,D3); 

% Each layer is treated separately:
for layer=1:D3
    % Transformeing the matrix to percentage:
    aa_freq = zeros (M,20);
    for m=1:M % m is the seq number
        aa_sum = sum(aa_total(m,:,layer),2);
        if aa_sum > 0
            % here we use elementwise division:
            aa_freq(m,:) = 100*aa_total(m,:,layer)./aa_sum;
        end
    end

    aa_freq_group = zeros (Ngroups,20);
    % To finde the mean frequence of the baseline group (1):
    % (mean in dim 1 specified in case of vector, ie. only one seq in the group)
    aa_freq_group(1,:)=mean(aa_freq(1:nGroup(1),:),1); 
    % To find the mean frequency composition in the other groups of sequences:
    startSeq=nGroup(1);
    for n=2:Ngroups
        stopSeq=startSeq+nGroup(n);
        aa_freq_group(n,:)=mean(aa_freq(startSeq+1:stopSeq,:),1); % dim 1 specified
        startSeq=stopSeq;
        % To returne data relative to baseline group:
        aaDfreq(n-1,:,layer) = aa_freq_group(n,:)-aa_freq_group(1,:); 
    end

end %layer loop
