% script to analyze TIGR double loop/double reference data

close all;
clear;

% load in data
load 'data_gpix.mat';

% put double loop and double reference together
pmtbig = [loopdata, refdata];

% make data object with replicates equal to 2
tigrdata = createData(pmtbig, 2);
tigrdata

% some diagnostic plot ...
% riplot
riplot(tigrdata);
close all;
% spatial pattern of the ratios
lratio = make_ratio(tigrdata);
arrayview(row, col, lratio);
close all;

% data normalization, use intensity lowess - spatial lowess is very slow
tigr_data_lowess = malowess(tigr.data, row, col, 'lowess');
save 'tigr_data_lowess' tigr_data_lowess;
%load tigr_data_lowess;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now start the data analysis - we will analyze the whole 
% data as one experiment 
% 
% STEP I: ANOVA and F test with:
%    - fixed effect model
%    - unrestricted residual shuffling
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% make a design object
design.sample = [loopvarid refvarid];
% note that sample 6 is reference sample
design.sample(design.sample==6) = 0;

% make model object
model_full_fix = makeModel(tigr_data_lowess, design, 'AG+DG+SG+sample');
model_full_fix
% ANOVA 
anova_full_fix = fitmaanova(tigr_data_lowess, model_full_fix);
anova_full_fix
% make null model
nullmodel = makeModel(tigr_data_lowess, design, 'AG+DG+SG');
% Do F test
ftest_fix = permTest(tigr_data_lowess, nullmodel, model_full_fix);
save 'ftest_fix' ftest_fix;
%load ftest_fix
ftest_fix

% volcano plot
idx_fix = volcano(anova_full_fix, ftest_fix);

% Clustering analysis
% bootstrap HC on samples
nodeobj = boothc(tigr_data_lowess, model_full_fix, idx_fix, 'sample', ...
  100, 'seuclidean', 'sample', 'vg');
% build consensus tree from bootstrap result
ct_sample_fix = consensus(nodeobj, 0.95);

% bootstrap Kmean on genes
[class grp] = bootkmean(tigr_data_lowess, model_full_fix, idx_fix, 'sample', ...
  5, 100, 0.7, 'gene', 'vg');
VGprofile(anova_full_fix.sample, grp);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP II: ANOVA and F test with:
%    - mixed effect model, AG and SG as random
%    - unrestricted residual shuffling
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model_full_mix = makeModel(tigr_data_lowess, design, '!AG+DG+!SG+sample');
model_full_mix

% ANOVA 
anova_full_mix = fitmaanova(tigr_data_lowess, model_full_mix);
save 'anova_full_mix' anova_full_mix;
%load anova_full_mix
anova_full_mix

% make null model and do f test
ftest_mix = test(tigr_data_lowess, model_full_mix, 'sample');
save 'ftest_mix' ftest_mix;
%load ftest_mix
ftest_mix

% calculate p values
ftest_mix = pval(ftest_mix, 4, 34);
ftest_mix
ftest_mix.F1

% volcano plot
figure;idx_mix = volcano(anova_full_mix, ftest_mix);