function [H, P]=propNormalityTest(q, aaOrder, alpha);
% propertyNormalityTest is used to find if frequences or substitutions 
% of a physicochemical quantity q has a Normal distribution or not by a 
% Kolmogorov-Smirnov type of goodnes of fit test

% Program by Steinar Thorvaldsen, steinart@math.uit.no, Okt. 2005. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last changes 31. Jan 2007.

% Use:   [Hnormal, Pvalue]=propNormalityTest(q, 2, 0.2);

% Input
% q:     Vector with 20 physicochemical properties, one for each amino acid.
% aaOrder: the order of comaprision (1=aa frequences, 2=aa substitutions)
% alpha: the level for the test. 
% Output
% H (normality): 0 (yes) or 1 (no)

if nargin < 3
    alpha=0.1; %default, may be in range 0.01-0.2 in the lillietest
end
if nargin < 2
    aaOrder=2; % standard substitutions
end

if aaOrder==1 %(aa frequences):
% The UniProt database distribution of amino acids:
Xpercent= [7.8 5.29 4.22 5.3 1.57 3.93 6.59 6.93 2.27 5.91 9.62 5.93 2.37 4.02 4.85 6.89 5.46 1.16 3.09 6.69];
% Generating a representative sample size for each amino acid
% with total sample size N ~50 (with roundoff effects):
aaSampleSize=round(Xpercent.*50./(sum(Xpercent))); %integer values       
N=sum(aaSampleSize);
% Then we sample N amino acids for final input in the test:
dsample=zeros(1,N);
% Sampling N (~100) data for the test:
n=0;
for i = 1:20
    if aaSampleSize(i) >= 1
        n2=n+aaSampleSize(i);
        dsample(n+1:n2)=q(i);
        n=n2;
    end
end
    
else %aaOrder=2 (substitutions):
% Mutation prob matrix (PAM1) scaled by 10^5:
% Ref.: DT Jones, WR Taylor and JM Thornton (1992). CABIOS vol 8:3 pp 275-282
PAM1= [ 0 27 24 42 12 23 66 129  5 19 28 22 11  6 99 264 267 1 4 193;
       41  0 19  8 21 125 20 102 74 13 34 390 10 3 36 69  38 18 8 11;
       43 23  0 284 6 31 36  58  92 26 12 150  8 3 6 344 137 0 23 11;
       63  8 235 0  2 21 478 95  24 6  6  17  4  1 6  40 25  1 15 21;
       44 52 13  5  0  4  3  41  17 8 15  3  10 28 6 147 28 16 68 41;
       43 154 33 27 2  0 211 17 130 4 64 176 11  2 81 37 31  2  8 12;
       82 16 25 398 1 140  0 83  6  6  9 103  4  2 10 21 19  2  2 31;
       135 70 33 66 11 10 70  0  5  3  6  16  3  2 11 129 19 8  2 32;
       17 164 171 53 15 233 15 15 0 10 49 31  8 18 58 51 28  2 189 8;
       28 12 21  6  3  3  7  4  4  0 212 12 113 31 5 28 149 2 10 630;
       24 19  6  3  3 29  6  5 12 122 0  9 90 101 53 40 16  8  8 117;
       28 334 108 14 1 122 107 20 12 11 13 0 15 1 11 32 57  1  3   8;
       36 22 14 10  8 19 11 10  8 253 350 37  0 18 8 19 123 3  6 201;
       11  3  3  2 14  2  3  4 11 41 230  1  10  0  8 65  8 8 179 40;
       150 36 5  7  3 66 12 16 26  5 97  13  4   6  0 190 69 1  4 14;
       297 51 214 30 44 22 19 139 17 21 54 28 7 38 140 0 278 4 20 27;
       351 33 100 22  9 21 20 24 11 134 25 57 49 6 59 325  0 1  6 76;
       7 65  1  3  23  7   7 41  3  7 49  5  5 22  4  21  5  0 24 16;
       11 12 30 23 43 10  4  4 134 16 22  5  4 222  6 43 12 11  0 11;
       226 9  7 16 13  7 29 35  3 504 161 7 71 24  11 28 67  3  5 0];
 
% Generating a representative sample size for each mutation
% with total sample size N ~100 (with roundoff effects):
MutSampleSize=round(PAM1.*136./(sum(sum(PAM1))));          
N=sum(sum(MutSampleSize));

% We also compute the change of q for each possible substitution:
dq=zeros(20,20);
for i=1:20
    for j=1:20
        dq(i,j)=q(i)-q(j);
    end
end

% Then we sample N differences (substitutions)
% for final input in the test:
dsample=zeros(1,N);
% Sampling N (=100) data for the test:
n=0;
for i = 1:20
    for j=1:20
        if MutSampleSize(i,j) >= 1
            n2=n+MutSampleSize(i,j);
            dsample(n+1:n2)=dq(i,j);
            n=n2;
        end
    end
end
end %if aaOrder

% lillietest: Lilliefors test for goodness of fit to a normal distribution:
[H,P] = lillietest(dsample,alpha); % Matlab statistical toolbox
% Alternative jbtest: Jarque-Bera test for goodness-of-fit to a normal distribution:
%[H,P] = jbtest(dsample,alpha); % Asyptotic test in Matlab stat toolbox
