function [prop_P,prop_name]=propPairedPvalues(Sgrouped, level, levelVar, threshPlotP);
% propertiesPairedPvalues is used to find and plot the significant changes
% between aligned sequences from the two extreme opposite groups in Sgrouped.
% The statistical test is data-adaptive based on a paired t-test or a Wilcoxon test.
% Remark: In this procedure the analyses is based on all sites in the alignment(s),
% and gaps are problematic. If you instead prefere to run a test based on the MEAN value of
% each sequence, you may apply the regression procedure propRegress.m or propRegressKendall.m
% Program by Steinar Thorvaldsen, steinart@math.uit.no, Nov. 2004. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last changes 26. July 2006.

% Use:     [prop_P,prop_name]=propPairedPvalues(Sgrouped, level, levelVar, 0.01);

% Input
% Sgrouped: Cell array containing the gene sequence from 2 or more
%           different groups of organisms. All seq must be of same length.
% level:    Vector containin information on 2D or 3D about the aligned sequences. 
%           The number of layers may be 1.
% levelVar: Variation level at each site in the sequence. Conserved sites has levelVar=1.
% threshPlotP: Threshold p-value for performing plots.
% A number,Nprop, of chemical property to filter are also read from a property table for aa.

% Output
% prop_P:    Nprop x layers array. Statistical P-values for the observed seq difference.
% prop_name: String array with name of properties

% Reconstructing the dimensions 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
[M1,N1] = size(cell2mat(Sgrouped(1))); %N1 must be same for all seq groups

% The function is made flexible conserning input parameters:
if nargin < 4
    threshPlotP=0; % default no plotting
end
if nargin < 3
    levelVar=2*ones(N1,1,'int8'); % no variation levels
end
if nargin < 2
    level=ones(N1,1,'int8'); % no structural levels
end

levelVar1=find(levelVar<=1); % the conserved sites
D3=max(level); % the number of layers in the level filter
if D3 > 1 % Constructe D3+1 cells with index vectors for each level:
    D3=D3+1; %append extra space for the sum of all other layers
    for n=1:D3-1
        level_index{n} = {find(level==n)};
    end
    level_index{D3} = {(1:N1)}; %containing all levels (should gaps be excluded?)
    %level_index{D3} = {find(level<=D3-2)}; %containing all template levels, gaps excluded
end 

% Read tabular values for 60+ amino acid properties from Excel file:
% NB: not all gap-values (no 21) may be NaN (then xlsread treats the column as text)
[prop_v,prop_name]=xlsread('aa_prop60.xls');
[Mp,Np]=size(prop_v);
prop_name=prop_name(3:Mp+2,22); % to extract names from column nr 22, when top margin=2
prop_P = zeros (Mp,D3); % the P-values for Mp properties separated in D3 layers

important_hit=[55]; %make a list after prescanning
for mp=1:Mp %important_hit %1:Mp %scanning through all properties
    % Perform a goodnes of fit test for normality of data or not:
    [Hnormal]=propNormalityTest(prop_v(mp,:), 0.2); %Hnormal=1(non-normal), Hnormal=0(normal)
    % Computing average property values along the two extreme seq groups:
    y_prop_upper = filterAlignedseq(cell2mat(Sgrouped(1)), prop_v(mp,:),1); % y-values along baseline seq group
    y_prop_lower = filterAlignedseq(cell2mat(Sgrouped(Ngroups)), prop_v(mp,:),1); %the oposite group
    diff_prop = y_prop_lower - y_prop_upper;
    
    % Each layer is analysed separately:
    for layer=1:D3
        current_ind=cell2mat(level_index{layer});
        current_ind_no_conserved=setdiff(current_ind,levelVar1); % fully conserved sites shold be excluded from t-test
        % The procedure is data-adaptive:
        if Hnormal ==0 %normal data
         [h,prop_P(mp,layer),ci] = ttest(diff_prop(current_ind_no_conserved)); % The paired t-test
        else %non-normal data requires a non-parametric test based on ranks:
           [prop_P(mp,layer),h] = signrank(diff_prop(current_ind_no_conserved)); % The Wilcoxon signed rank test
        end
        
        if layer ~= D3-1 & prop_P(mp,layer) < threshPlotP %0.01 % candidate for an interesing result
            figure
            stem (diff_prop(current_ind));
            title(['Difference in prop ', num2str(mp),' between populations, level ', num2str(layer),', p=',num2str(prop_P(mp,layer))]);
            hold off
            
            %We also make a smoothed plot with all groups and a broader filter window (w=3):
            figure            
            y_prop_level = zeros(Ngroups,length(current_ind));
            for ng=1:Ngroups
                Sg=cell2mat(Sgrouped(ng));
                y_prop_level(ng,:) = filterAlignedseq(Sg(:,current_ind), prop_v(mp,:),3); %amino acids at current level
                plot (y_prop_level(ng,:));
                hold all
            end
            title(['Property ', num2str(mp),' variation in populations (w=3) at level ', num2str(layer),', p=',num2str(prop_P(mp,layer))]);
            hold off         
        end
        if nanmean(diff_prop(current_ind_no_conserved)) < 0 % nanmean skips all NaN-values
            prop_P(mp,layer) = -prop_P(mp,layer); % just to note decreasing trends by a minus
        end
    end
        
end
