% the script to analyze Bev Paigen's 28 array experiment
% Only 300 genes are included in this example for the sake of speed
% The data are after-normalization. The text files are output from
% R/maanova

clear;
close all;

% read in data 
[data, varnames, geneid] = tblread('paigen300.txt', '\t');
% get pmt data
pmt = data(:,5:end);

% read in design information
[design.Strain, design.Diet, design.sample] = textread('paigen300design.txt', ...
  '%s%s%d', 'headerlines', 1);

% make data object
paigen300data = createData(pmt, 2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fixed model analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make several model objects and fit ANOVA model
% make full model object
model_full_fix = makeModel(paigen300data, design, 'DG+AG+SG+Strain+Diet+Strain*Diet');
anova_full_fix = fitmaanova(paigen300data, model_full_fix);
% make model object without interaction
model_noint_fix = makeModel(paigen300data, design, 'DG+AG+SG+Strain+Diet');
anova_noint_fix = fitmaanova(paigen300data, model_noint_fix);
% make model object without strain effect
model_nostrain_fix = makeModel(paigen300data, design, 'DG+AG+SG+Diet');
anova_nostrain_fix = fitmaanova(paigen300data, model_nostrain_fix);
% make model object without diet
model_nodiet_fix = makeModel(paigen300data, design, 'DG+AG+SG+Strain');
anova_nodiet_fix = fitmaanova(paigen300data, model_nodiet_fix);

% permutation test - all tests use non-restricted residual shuffling
% test for interaction effect
test_int_fix = permTest(paigen300data, model_noint_fix, model_full_fix);
save 'test_int_fix' test_int_fix;
%load test_int_fix;
idx_int_fix = volcano(anova_full_fix, test_int_fix);
title('Volcano plot for interaction test - fixed model');
% test for strain effect 
test_strain_fix = permTest(paigen300data, model_nostrain_fix, model_noint_fix);
save 'test_strain_fix' test_strain_fix;
%load test_strain_fix;
figure; idx_strain_fix = volcano(anova_noint_fix, test_strain_fix);
title('Volcano Plot for Strain effect test - fixed model');
% test for strain effect 
test_diet_fix = permTest(paigen300data, model_nodiet_fix, model_noint_fix);
save 'test_diet_fix' test_diet_fix;
%load test_diet_fix;
figure;idx_diet_fix = volcano(anova_noint_fix, test_diet_fix);
title('Volcano Plot for Diet effect test - fixed model');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% mixed model analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make several model objects and fit ANOVA model
% make full model object
model_full_mix = makeModel(paigen300data, design, 'DG+!AG+!SG+Strain+Diet+Strain*Diet+!sample');
anova_full_mix = fitmaanova(paigen300data, model_full_mix);
save 'anova_full_mix' anova_full_mix;
%load anova_full_mix;

% model object without interaction
model_noint_mix = makeModel(paigen300data, design, 'DG+!AG+!SG+Strain+Diet+!sample');
anova_noint_mix = fitmaanova(paigen300data, model_noint_mix);
save 'anova_noint_mix' anova_noint_mix;
%load anova_noint_mix;

% test interaction effect
test_int_mix = test(paigen300data, model_full_mix, 'Strain*Diet');
save 'test_int_mix' test_int_mix;
%load test_int_mix;
% calculate the P valus, numerator's df is 2, denominator's df is 6
test_int = pval(test_int_mix, 2, 6);
% volcano plot
figure; idx_int_mix = volcano(anova_full_mix, test_int);
title('Volcano Plot for Interaction effect test - mixed model');

% Use model without interaction to test Strain and Diet effects
% Strain effect
test_strain_mix = test(paigen300data, model_noint_mix, 'Strain');
save 'test_strain_mix' test_strain_mix;
%load test_strain_mix;
test_strain = pval(test_strain_mix, 2, 8);
figure; idx_strain_mix = volcano(anova_noint_mix, test_strain);
title('Volcano Plot for Strain effect test - mixed model');

% Diet effect
test_diet_mix = test(paigen300data, model_noint_mix, 'Diet');
save 'test_diet_mix' test_diet_mix;
%load test_diet_mix;
test_diet = pval(test_diet_mix, 1, 8);
figure; idx_diet_mix = volcano(anova_noint_mix, test_diet);
title('Volcano Plot for Diet effect test - mixed model');