% DeltaProt1strata: a protein comparision toolbox, main program sample script. 
% A comparative analysis of ONE aligned protein dataset with 3 taxionomic strata
% Program by Steinar Thorvaldsen, steinar.thorvaldsen@uit.no, Nov. 2004. 
% Ref.: http://www.math.uit.no/bi/deltaprot/
% S.Thorvaldsen and E.Ytterstad: "Environmental adaptation of proteins: 
% Regression models with simple physicochemical properties".
% Computational Biology and Chemistry. Vol. 33 (5) 2009, Pages 351-356.

% Last uppdated 22.Oct. 2010

% Input:
% One file with m aligned protein sequences extracted from different organisms.
% (fasta format). The file must be stored in the subdirectory \data1\
% To each organisme a biological or phenotype parameters may be attached 
% to the header of the sequence:
% Topt[C]:  Optimal growth temperature for the organisme
% Taxionomic strata:   Order of organisems included in the model using dummy varables

% Auxilary predictions to be applied beforehand (external servers)
% with results stored at the top or bottom of the alignment:
% GETAREA:  (http://www.scsb.utmb.edu/cgi-bin/get_a_form.tcl) to find the
% Water Accessible Surface Area.
% (0-completely buried, 9-fully exposed) 
% If none of the sequences has solved structure in 3D, the program SABLE
% may predicte: (http://sable.cchmc.org/)

% Procedures in DeltaProt1 Comparision:

% I-III: properties tests
%   Comparing physiochemical properties in sequence pairs and groups.
%   We apply 3 statistical tests (Ordinary regression, Wilcoxon paired, one-way ANOVA).

% IV: amino acid composition
%   Input: protein data set, and optional vector with structural level 2D and 3D.
%   Output matrix: groups x 20 x levels. Amino acid frequency matrix.

clear all;
%---------------------- Read and convert data --------------------------
aa_seq = dir('./data1/*.fas'); % reads the filname from data directory into a structure array with fields
fileNo=size(aa_seq,1); % the number of files in the directory
if fileNo ~= 1
    error ('The number of files in directory /data1 is different from 1. Use script DeltaProt.m for multiple sequence analysis');
end
top_lines=4; %line 1..4 in alignment reserved for 2D/3D information
% L is the number of lines in the alignment file, n is the length, 
% S is the character sequence array, and names are the headers:
[L, n, S, names] = readfasta(strcat('./data1/',aa_seq(fileNo).name)); %strcat adds the two strings
Sint=aa2int(S(top_lines+1:L,:)); % integer number representation of alignment S
siteVar = findSeqSiteVar(Sint); % finding the variation and conserved sites (Var=1) along ALL sequences
% Computational filters added at the top or bottom of the alignment:
level2D=dssp2int(S(1,:)); % secondary structure(Alpha; Beta Loop)
level3D=wasa2int(S(3,:)); % internal/external amino acids(0..9)
%levelDNA=DNAdist2int(S(35,:)); % closeness to DNA (3=8Å, 2+3=10Å, 1+2+3=12Å)

%------------I. Inferring group trends based on properties---------------
% For stratified regression: Vibrionales used as baseline + two dummy variables:
X = ones (L-top_lines,3)*NaN; % the regressor variable
X(:,1)=[1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3 2 3]'; %temp habitat 1-3
X(:,2)=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0]'; %Alteromonadales
X(:,3)=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1]'; %Desulfobacterales

% We look for significant physiochemical trend changes in the groups:
[b,prop_P,prop_name]=propRegress(Sint, X, level3D, siteVar);
 % FDR correction for multiple testing:
yer =[4 5 7 8 30 31 36 47 54 55 56 57 58 59 85]; %selected physicochemical properties
xer = [1 2 3 5]; %avoiding level 4
FDR (abs(prop_P(yer,xer)), 0.05)

%-----------II. Compute and test sequence properties by paired tests -----
% The three temperature habitats:
t31_40=[1 2 3 4 5 6 17 18 19]; %mesophiles
t21_30=[7 8 9 10 11 12 20 21 22 23 24 25 34]; %intermediate
t1_20=[13 14 15 16 26 27 28 29 30 31 32 33 35]; %psycrophile
% We use cell arrays as a storage mechanism for dissimilar size of data:
Sgrouped= {Sint(t31_40,:); Sint(t21_30,:); Sint(t1_20,:)};
[paired_p2D,prop_name]=propPairedPvalues(Sgrouped, level2D, siteVar, 0.001);
%[paired_p3D,prop_name]=propPairedPvalues(Sgrouped, level3D, siteVar, 0.001);

%-----------III. Compute and test sequence properties by one-way ANOVA --
%[ANOVA_p2D,prop_name]=propAnovaPvalues(Sgrouped, level2D, siteVar);
%[ANOVA_p3D,prop_name]=propAnovaPvalues(Sgrouped, level3D, siteVar);

%--------------------- IV. Compute and plot aa-composition----------------
% Secondary structure:
%[aaComp2D,aaStd2D, aaTotal2D] = aaFreq(Sgrouped, level2D);
%aaFreqPlot(aaComp2D, aaStd2D, 1);
% Structural 3D location:
[aaComp3D,aaStd3D,aaTotal3D] = aaFreq(Sgrouped, level3D);
aaFreqPlot(aaComp3D, aaStd3D, 3);
