% DeltaProt: multiple protein comparison, main program sample script. 
% A comparative analysis of aligned, homolog protein sequences.
% Program by Steinar Thorvaldsen, steinart@math.uit.no, Nov. 2004.
% Ref.: http://www.math.uit.no/bi/deltaprot/
% Last uppdated 22.Oct. 2010

% Input:
% Fasta files with m aligned protein sequences extracted from different organisms.
% The files must be stored in fasta format in a subdirectory like \data\
% To each sequence some phenotypical parameters may be attached 
% to the header of the sequence:
% Topt[C]:  Optimal growth temperature for the organisme
% GC[%]:    The amount og G+C nucleotides in the genome
% pHopt:    Optimal pH for the organisme
% Popt[atm]:Optimal pressure for the growth of the organisme
% NaClopt[mM]:Optimal salt consentration for growth of the organisme
% Phylum:   Life Domain
% .... and so on.
% This parameter is used to group the data in an ordinal variable.
% NB: The grouping in ordinal populations is done manually in the script below.

% Auxilary predictions to be used 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)
% If none of the sequences has solved structure in 3D, another program may be used:
% SABLE:    Relative solvent accessibility of the residues
%           (http://sable.cchmc.org/)
%           (0-completely buried, 9-fully exposed) 
%           Confidence level of the RSA prediction (0-low confidence, 9-highly confident) 
% DSSP:     Finding secondary structure of sequences from structure
% If none of the sequences has solved structure in 3D, another program may be used:
% PsiPred:  Predicting secondary structure of sequence (H, E, C code) 
%           (http://bioinf.cs.ucl.ac.uk/psipred/psiform.html)
%           Confidence level of the SS prediction (0-low confidence, 9-highly confident) 

% Procedures in DeltaProt Comparision:
% I. aaDeltaFreq
%   Input: protein data set, and optional vector with structural level 2D and 3D.
%   Output matrix: groups x 20 x levels. Amino acid frequency matrix.

% II. substPairsCount
%   Input aligned protein sequence pair
%   Ouput substitution matrix with optional gaps (insertions/deletions)

% IIIa. propAnovaPvalues
%   Comparing sequence properties by one-way ANOVA
% IIIb. propPairedPvalues
%   Comparing physiochemical properties in sequence groups.
%   Visualise data by using: filterAlignedseq
% IIIc. propRegressKendall
%   Comparing poulation properties by nonparametric regression

clear all;
warning off
% You may use the global constant 'maxlevel' to force the analyses of all
% genes to be done with the same number of 2D and 3D levels:
% (not relevant when only one alignment is to be analysed)
global maxlevel;
    maxlevel=4; % levels of 2D/3D in analyses, may be changed ****
    
%---------------------- Read and convert data --------------------------
aa_seq = dir('./data/*.stn'); % reads the filname from data directory into a structure array with fields
maxNo=size(aa_seq,1); % the number of files in this particular case
M=12; %the number of data lines represented in this data set (the fasta-files)

% Defining sequence populations of all 7 sequences (grouped by temperature):
maxPop=3;
t31_40=[2 3 4]; % the mesophile sequences, nr 1 is a copy
t21_30=[5];
t1_20=[6 7]; % the psychrophile sequences.
extra_lines=5; %line 8..12 reserved for 2D/3D information and optional template seqence

%Here dim 2 and 4 must match with data! May be changed:
aaDeltafreq2D = zeros (maxNo,maxPop-1,20,maxlevel+1); % for storing changes in frequencies
aaDeltafreq3D = zeros (maxNo,maxPop-1,20,maxlevel+1); % for storing changes in frequencies
proteinBlock = zeros(maxNo,1); % for storing the length of each protein sequence block
Stotal=[];
siteVartotal=[];
level2Dtotal=[];
level3Dtotal=[];

for fileNo=1:maxNo
    % 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('./data/',aa_seq(fileNo).name)); %strcat adds the two strings
    proteinBlock(fileNo) = n; % length of this alignment
    Sint=aa2int(S(1:L-extra_lines,:)); % integer number representation of alignment S
    siteVar = findSeqSiteVar(Sint); % finding the variation and conserved sites (Var=1) along ALL sequences
    Sgrouped= {Sint(t31_40,:); Sint(t21_30,:); Sint(t1_20,:)};
    % Computational filters added at the top or bottom of the alignment:
    level2D=dssp2int(S(9,:)); % secondary structure(1,2,3)
    level3D=wasa2int(S(11,:)); % internal/external amino acids(only levels 1,2,3)
    Stotal=[Stotal,Sint]; % store data by just extending the arrays in dim 2
    siteVartotal=[siteVartotal,siteVar];
    level2Dtotal=[level2Dtotal,level2D]; 
    level3Dtotal=[level3Dtotal,level3D];
    aaDeltafreq2D(fileNo,:,:,:) = aaDeltaFreq(Sgrouped, level2D);
    aaDeltafreq3D(fileNo,:,:,:) = aaDeltaFreq(Sgrouped, level3D);
    
end
%---------------------- I. Plot the frequence variation of amino acids---------
% Secondary structure:
%aaDeltaFreqPlot(aaDeltafreq2D);
% Structural 3D location:
aaDeltaFreqPlot(aaDeltafreq3D);


%---------------------- II. Compute and plot aa-substitution matrix-------
% Secondary structure of protein (exsample)
substMP2D = substPairsCount(Stotal(3,:), Stotal(6,:), level2Dtotal, 1);
PvalueMPvsPM2D = substPairsPvalues(substMP2D);
%substPairsPlot(substMP2D,PvalueMPvsPM2D, 0.05);
substPvalueSort(substMP2D,PvalueMPvsPM2D);
substCountSort(substMP2D,PvalueMPvsPM2D);

% 3D structure of protein (exsample)
substMP3D = substPairsCount(Stotal(t31_40,:), Stotal(t1_20,:), level3Dtotal, 0, 1);
PvalueMPvsPM3D = substPairsPvalues(substMP3D);
substPairsPlot(substMP3D,PvalueMPvsPM3D, 0.01);
substPvalueSort(substMP3D,PvalueMPvsPM3D);
substCountSort(substMP3D,PvalueMPvsPM3D);
substMM3D = substPairsCount(Stotal(t31_40,:), Stotal(t31_40,:), level3Dtotal, 0, 0);
PvalueMPvsMM3D = substPairsPvalues(substMP3D,substMM3D);


%-----------IIIa. Compute and test sequence properties by one-way ANOVA --
Sgrouped= {Stotal(t31_40,:); Stotal(t21_30,:); Stotal(t1_20,:)};
%[ANOVA_p2D,prop_name]=propAnovaPvalues(Sgrouped, level2Dtotal, siteVartotal);
%[ANOVA_p3D,prop_name]=propAnovaPvalues(Sgrouped, level3Dtotal, siteVartotal);

%-----------IIIb. Compute and test population properties by paired tests ----
% As in the first case we use cell arrays to store dissimilar size of data:
Sgrouped= {Stotal(t31_40,:); Stotal(t21_30,:); Stotal(t1_20,:)};
% We look for significant physiochemical differences in the two extreme groups:
%[paired_p2D,prop_name]=propPairedPvalues(Sgrouped, level2Dtotal, siteVartotal);
[paired_p3D,prop_name]=propPairedPvalues(Sgrouped, level3Dtotal, siteVartotal);

%-----------IIIc. Compute and test poulation properties nonparametric regression---------
Sgrouped= {Stotal(t31_40,:); Stotal(t21_30,:); Stotal(t1_20,:)};
% We look for significant physiochemical trend changes:
%[b2D,cumKendall_p2D,prop_name]=propRegressKendall(Sgrouped, level2Dtotal, siteVartotal, proteinBlock);
[b3D,cumKendall_p3D,prop_name]=propRegressKendall(Sgrouped, level3Dtotal, siteVartotal, proteinBlock);
% For ploting slopes of one property (No 5):
Property5_Surfacelevel3_Allprot_Slope=sort(squeeze(b3D(5,3,:,2))); 
figure;
bar (Property5_Surfacelevel3_Allprot_Slope);
warning on
