function [p, t,df] = mattest(X1, X2, varargin)
%MATTEST Two-sample t-test for identifying differentially expressed genes
% from two different conditions/phenotypes.
%
%   P = MATTEST(X1,X2) performs a t-test for differential expression with a
%   standard two-sided t-test on every row (gene) of the expression
%   matrices, and returns the results of the p-values for every gene. The
%   data are assumed to be normally distributed in each class with equal
%   variances. X1 and X2 must have the same number of rows. 
%
%   [P, T] = MATTEST(...) returns the t-scores T for each gene.
%
%   [P, T, DF] = MATTEST(...) returns the degrees of freedom DF of the
%   test.
%
%   MATTEST(...,'PERMUTE', P) perform permutation tests to if P is set to
%   TRUE.  The default number of permutation is 1000. P can also be the
%   number of permutations if is set to be greater than 2. The default is
%   FALSE.
% 
%   MATTEST(...,'BOOTSTRAP', B) perform bootstrap tests to if B is set to
%   TRUE.  The default number of bootstrap is 1000. B can also be the
%   number of bootstrap if is set to be greater than 2. The default is
%   FALSE.
% 
%   MATTEST(...,'SHOWHIST',TRUE) displays histograms of the p-value
%   distribution and t-score distribution. 
%
%   MATTEST(...,'SHOWPLOT',TRUE) displays normal t-scores quantile plot.
%   The points with t-score greater than the (1 - 1/(2N)) quantile or less
%   than the 1/(2N) quantile are marked with red circles. N is total number
%   of genes in the t-test.
% 
%   MATTEST(...,'LABELS',LABELS) allows you to specify a cell array of gene
%   labels or names for each row of the input data. If LABELS are defined,
%   then clicking on a data point on the plot will show the LABEL
%   corresponding to that data point.
% 
%   Example:
%       load prostatecancerexpdata; 
%       [p,t]=mattest(dependentData,independentData, 'showplot', true);
% 
%       % Permutation t-tests.
%       [p,t]=mattest(dependentData,independentData, 'permute', true,...
%                           'Showhist', true', 'showplot', true);
% 
%       % Bootstrap t-tests.
%       [p,t]=mattest(dependentData,independentData, 'bootstrap', 2000,...
%                           'Showhist', true', 'showplot', true);
% 
%   See also GCRMA, MAFDR, MAIRPLOT, MALOGLOG, MAPCAPLOT, MAVOLCANOPLOT,
%   RMASUMMARY.

% Copyright 2005-2007 The MathWorks, Inc.
% $Revision: 1.1.10.7 $   $Date: 2007/11/22 16:39:12 $ 

% References: 
%  [1] W. Huber, A. von Heydebreck, H. Sültmann, A. Poustka, M. Vingron;
%      "Variance stabilization applied to microarray data calibration and
%      to the quantification of differential expression", Bioinformatics
%      (2002) 18 Suppl1 S96-S104.


bioinfochecknargin(nargin,2,mfilename)

%== Parse input options
appdata = parse_inputs(varargin{:});

appdata = check_inputdata(X1, X2, appdata);

% do ttest
if appdata.nperms ~= 0
    [p, t, df] = permutationtests(X1, X2, appdata.nperms);
elseif appdata.nboots ~= 0
    [p, t, df] = bootstraptests(X1, X2, appdata.nboots);
else
    [h_dum,p,ci_dum,stats] = ttest2(X1, X2, [],[],[],2); %#ok
    t = stats.tstat;
    df = stats.df;
end

if appdata.qplotFlag  
    hqpfig = plotquantilet(t, df, appdata.labels); %#ok
end

if appdata.histFlag
    appdata.hdfig = plothistogram(t, p); %#ok
end
end

%*********** Callbacks and helper functions ***************
function [p, t, df] = permutationtests(X1, X2, nperms)
% Perform permutation tests
[p, t, df] = resampling(X1, X2, nperms, false);
end

%------------------------------------------------------
function [p, t, df] = bootstraptests(X1, X2, nboots)
% Perform bootstrap two-sample ttest. 
%   Resampling under H0 
%   Mix the two samples
%   Draw bootstrap sample from each group
%   Compute p-value = sum(I(t_b >= t_0))/B.

[p, t, df] = resampling(X1, X2, nboots, true);
end

function [p, t, df] = resampling(X1, X2, nsamples, replace)
% Number of samples
n1 = size(X1, 2);
n2 = size(X2, 2);
n = n1 + n2;
X = [X1 X2];
m = size(X1, 1);

% Compute observed t statitics
t = computeGroupT(X1, X2);

% Conduct a permutation (bootstrap) test to create a sample from the null
% distribution of the t statistics. Permute the transcripts at random (with
% replacement for bootstrap), compute t-scores among the random partitions
% of the actual samples, repeat for a large number permutations.
% Bin the samples
edge = max(findedge(t), 20);
t0edges = linspace(-edge, edge, getUpperGeneNumber(m)+1)'; % #bins odd -> edge at zero
t0bins = t0edges(1:end-1) + diff(t0edges)/2;
t0 =zeros(size(t0edges));

for i = 1:nsamples
    if replace % with replacement for Bootstrap
        grp1 = X(:, randsample(n, n1, true));
        grp2 = X(:, randsample(n, n2, true));
    else % without replacement for permutation
        [grpidx1, grpidx2] = crossvalind('LeaveMOut', n, n2);
        grp1 = X(:, grpidx1);
        grp2 = X(:, grpidx2);
    end
    t0rep = computeGroupT(grp1, grp2);
    t0 = t0 + histc(t0rep, t0edges);
end

t0 = t0(1:end-1); % drop the degenetated uppermost bin

% Estimate the null distribution of |t| and compute its empirical
% cumulative distribution function.
[f0, x0] = ecdf(abs(t0bins), 'freq', t0);
x0(1) = 0;

% Estimate p-values corresponding to truly null hypotheses for the observed
% t-scores from the t-test
pp=interp1(x0, f0, 'pchip', 'pp');
p=abs(1-ppval(pp, abs(t)));

%== Degree of freedom
df = computeDF(n1, n2, t);
end

%--------------------
function [t, grpmean, grpvarm] = computeGroupT(X1, X2)
% Return group means var means and t scores
ngrp1 = size(X1, 2);
ngrp2 = size(X2, 2);
grpmean = [mean(X1,2) mean(X2,2)];
grpvarm = [var(X1,[],2)./ngrp1 var(X2,[],2)./ngrp2];
t = -diff(grpmean,[],2) ./ sqrt(sum(grpvarm,2));
end

%----------------------------------------------------------
function df = computeDF(n1, n2, t)
% Reture degree of freedom
df = n1 + n2 -2;
if ~isscalar(t)
    df = repmat(df,size(t));
end
end

%----------------------------------------------------------
function hqpfig = plotquantilet(tscore, df, labels)
% Do quantile plot
% df - degree of freedom.

hqpfig = figure('Units', 'Normalized',...
                'Tag', 'ttestquantileplot',...
                'Visible', 'off');

% Remove the NaNs 
goodidx = ~isnan(tscore);
tscore = tscore(goodidx);
df = df(goodidx);

if ~isempty(labels)
    labels = labels(goodidx);
end

[tscore, ord] = sort(tscore);
if ~isempty(labels)
    labels = labels(ord);
    
    if ~ischar(labels) && isnumeric(labels)
       labels = cellstr(num2str(labels(:)));
    end
    
end

% number of genes
N = numel(tscore);
% Get quatile function of t stats
Xq = (1:N)';
Xq = (Xq - 0.5)/N;
tquantile = tinv(Xq, df);

haxes = axes;
set(haxes, 'Parent', hqpfig,...
           'Box', 'on',...
           'GridLineStyle', '--');

if tquantile ~= 0
    set(haxes, 'XLimMode', 'manual',...
               'YLimMode', 'auto',...
               'XLim',[min(tquantile), max(tquantile)]);
end

%up-regulated threshold
upThd = tinv( 1- 0.5/N, df);
upIdx = find(tscore > upThd);
%down-regulated threshold
dnThd = tinv( 0.5/N, df);
dnIdx = find(tscore < dnThd);

hplot = plot(tquantile, tscore, 'b.',...
             tquantile(upIdx), tscore(upIdx),'ro',...
             tquantile(dnIdx), tscore(dnIdx), 'ro',...
             tquantile, tquantile, 'k-.');

set(hplot, 'parent', haxes);
set(hplot(1), 'Tag', 'QuantileLine',...
              'DisplayName', 'Quantile');
if ~isempty(upIdx)
  set(hplot(2), 'Tag', 'UpQuantile',...
                'DisplayName', 'Significant');
end

if ~isempty(dnIdx)
    set(hplot(end-1), 'Tag', 'DownQuantile',...
                      'DisplayName', 'Significant')
end

set(hplot(end), 'DisplayName', 'Diagonal')


%== Turn off data cursor for other lines
hb=hggetbehavior(hplot(end), 'Datacursor');
set(hb, 'Enable', false)

%== Enable only the quantile line
dcm_obj = datacursormode(hqpfig);
set(dcm_obj,'UpdateFcn',{@locaDataCursorUpdate, dcm_obj, labels, upIdx, dnIdx})

% Set plot labels
xlabel('Theoretical quantile')
ylabel('Sample quantile')
title('Normal Quantile Plot of t')
set(hqpfig, 'Visible', 'on');
end

function txt = locaDataCursorUpdate(hobj,eobj, dcm, labels, upIdx, dnIdx) %#ok
% Update the datacursor text.
tg = get(eobj, 'Target');
c=getCursorInfo(dcm);
if strcmpi(get(tg, 'Type'), 'line')
    type = get(tg, 'Tag');
    pos = get(eobj,'Position');
    if strcmpi(type, 'QuantileLine')            
        txt = {['Label: ', labels{c.DataIndex}],...
            ['T. quantile: ',num2str(pos(1))],...
            ['S. quantile: ',num2str(pos(2))]};
    elseif strcmpi(type, 'UpQuantile')
         txt = {['Label: ', labels{upIdx(c.DataIndex)}],...
            ['T. quantile: ',num2str(pos(1))],...
            ['S. quantile: ',num2str(pos(2))]};
    elseif strcmpi(type, 'DownQuantile')
         txt = {['Label: ', labels{dnIdx(c.DataIndex)}],...
            ['T. quantile: ',num2str(pos(1))],...
            ['S. quantile: ',num2str(pos(2))]};
    end
else
    txt = [];
    return;
end
end

%**************************************
function hdfig = plothistogram(tscore, pvalue)
% Plot histograms of tscores and p-values
hdfig = figure('Units', 'Normalized',...
                'Tag', 'ttesthistogram',...
                'Visible', 'off');
subplot(1,2,1);
bw = .5;
edge = findedge(tscore);
edges = -edge:bw:edge;
hgt1 = histc(tscore,edges);
bar(edges,hgt1,'histc');
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[1.0 0.3 0.3], 'EdgeColor', 'w')
xlim([-11 11]);
title('t-scores')
xlabel('t-score');
ylabel('Frequency')

subplot(1,2,2)
hist(pvalue, 0.025:0.05:.975);
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.5 0.5 0.8], 'EdgeColor', 'w')
title('p-values')
xlabel('p-value');
ylabel('Frequency')

suptitle('Histograms of t-test Results')
set(hdfig, 'visible', 'on')
end

%-----------------------------------------------
function n = getUpperGeneNumber(n)
% Return the rounded upper gene number
% for example, 3893 genes, return 4000

fn = floor(log10(n));
cn = ceil(log10(n));
if cn-fn == 0
    return;
end

m = ceil(n/(10^fn));
n = m*10^fn;
end
%---------------------------------------------------
function edge = findedge(t)
% return the rounded max |t|
maxt = getUpperGeneNumber(abs(max(t)));
mint = getUpperGeneNumber(abs(min(t)));
edge = max(maxt, mint);
end

%--------------------------------
function appdata = check_inputdata(X1, X2, appdata)
% Check input data type is numerical and contain the same number of genes
% or rows
if ~isnumeric(X1) || ~isreal(X1) || ~isnumeric(X2) || ~isreal(X2)
   error('Bioinfo:mattest:ExpressionValuesNotNumericAndReal',...
         'Expression values must be numeric and real.') 
end

% validate required inputs 
if size(X1,1)~=size(X2,1)
    error('Bioinfo:mattest:differentSize',...
        'X1 and X2 must have same number of rows.')
end

% Check the labels size
if(~isempty(appdata.labels))
    if ~isvector(appdata.labels) || numel(appdata.labels) < size(X1,1)
        warning('Bioinfo:mattest:MismatchedLabels',...
            ['The number of labels did not match the number of rows in input data.',...
            '\nThe quantile plot of t was not able to show the datatip.']);
        appdata.labels = [];
    end
    
    appdata.labels = appdata.labels(:);
else
    appdata.labels = (1:size(X1,1))';
end

% Check permutation and bootstrap options
if appdata.nboots ~= 0 && appdata.nperms ~= 0
     warning('Bioinfo:mattest:UsingPermuteOptionOnly',...
            'Using PERMUTE option only. BOOTSTRAP option will not be used.');
     appdata.nboots = 0;
end
end
%--------------------------------
function inputStruct = parse_inputs(varargin)
% Parse input PV pairs.

% Check for the right number of inputs
if rem(nargin,2)== 1
    error('Bioinfo:mattest:IncorrectNumberOfArguments', ...
        'Incorrect number of arguments to %s.', mfilename)
end

% Allowed inputs
okargs = {'permute', 'bootstrap', 'showhist', 'showplot', 'labels'};

% Defaults
inputStruct.qplotFlag = false;  % quantile plot
inputStruct.histFlag = false;   % histgram plot
inputStruct.labels = [];        % Gene labels
inputStruct.nperms = 0;         % number of permutation
inputStruct.nboots = 0;         % number of bootstrap

for j=1:2:nargin
    [k, pval] = pvpair(varargin{j}, varargin{j+1}, okargs, mfilename);
    switch(k)
        case 1 % permutation
            if opttf(pval)
                if isnumeric(pval)
                    if isscalar(pval)
                        inputStruct.nperms = ceil(double(pval));
                        if inputStruct.nperms <= 2
                            error('Bioinfo:mattest:InvalidPermutationNumber',...
                                'Number of permutation tests must be greater than 2.');
                        end
                    else
                        inputStruct.nperms = 1000;
                        warning('Bioinfo:mattest:PermuteNumberNoScalar',...
                            ['PERMUTE must be a single number or a logical.'...
                            '\nRun permutation test 1000 times by default.'])
                    end
                else
                    inputStruct.nperms = 1000;
                end
            else
                inputStruct.nperms = 0;
            end
        case 2 % bootstrap
            if opttf(pval)
                if isnumeric(pval)
                    if isscalar(pval)
                        inputStruct.nboots = ceil(double(pval));
                        if inputStruct.nboots <= 2
                            error('Bioinfo:mattest:InvalidBootstrapNumber',...
                                'Number of bootstrap tests must be greater than 2.');
                        end
                    else
                        inputStruct.nboots = 1000;
                        warning('Bioinfo:mattest:BootstrapNumberNoScalar',...
                            ['BOOTSTRAP must be a single number or a logical.',...
                            '\nRun bootstrap test 1000 times by default.'])
                    end
                else
                    inputStruct.nboots = 1000;
                end
            else
                inputStruct.nboots = 0;
            end
        case 3 % showhist
            inputStruct.histFlag = opttf(pval,okargs{k},mfilename);
        case 4 % show quantile plot
            inputStruct.qplotFlag = opttf(pval,okargs{k},mfilename);
        case 5 % labels
            inputStruct.labels = pval;
    end
end
end
