% Function: This program takes the consistency alpha cutoff and corrects
% for empirically measured technical error measured between the replicates

% Inputs:
% RepPValueCell - this is a cell array, where each entry contains the cell 
%    array with the replicate comparison p-values
% pCutoff - the p-value that is being corrected
% NRepValues - the number of replicates across the entire RepPValueCell,
%    entered as a numerical entry
% Filter - either "Y" or "N", to determine whether bad well were eliminated
%    from the false positive rate correction process because of obvious 
%    technical errors, defined when the median error in a well is less than
%    MHpvalue

% Outputs:
% alphaCutoff - the empirically adjusted p-value, based on the replicate
%    FPR
% RepPList4Fit - the sorted list of the replicate p-values used to make the 
%    correction for the alpha cutoff
% BadProbeCell - the list of wells that have obvious technical error

function [alphaCutoff, RepPList4Fit,BadProbeCell] = ...
    EmpiricalAlpha(RepPValueCell,pCutoff,NRepValues,Filter)

RepPList = zeros(1,NRepValues);
RepPListFiltered = [];
BadProbeCell = {};
ExpNum = length(RepPValueCell);
Index_Start = 1;

for i = 1:ExpNum

    CondNum = length(RepPValueCell{i});
    
    for j = 1:CondNum
        
        Index_End = Index_Start + length(RepPValueCell{i}{j}(:))-1;
        
        RepPList(Index_Start:Index_End) = RepPValueCell{i}{j};
        
        [Probe_N,IP_N] = size(RepPValueCell{i}{j}); 
        
        for alpha = 1:Probe_N

            Probea_MedError = median(RepPValueCell{i}{j}(alpha,:));

            if  Probea_MedError > 0.05/IP_N
                RepPListFiltered = [RepPListFiltered,RepPValueCell{i}{j}(alpha,:)];
            
            else
                BadProbeCell{end+1} = strcat('Exp',num2str(i),'_Cond',...
                    num2str(j),'_Probe#',num2str(alpha),'_MedianError=',...
                    num2str(Probea_MedError));
     
            end
            
        end
        
        Index_Start = Index_End+1;
        
    end
    
end

Sort_RepPList = sort(RepPList);

Sort_RepPList(isnan(Sort_RepPList)) = [];

if Filter == 'N'
    RepPList4Fit = Sort_RepPList;
else
    RepPList4Fit = sort(RepPListFiltered);
end
    
Sort_Index = round(pCutoff.*length(RepPList4Fit));

SIndex_m2 = min(max(Sort_Index-2,1),length(RepPList4Fit)-4);
SIndex_p2 = SIndex_m2+4;

PredictedFPR = (1:length(RepPList4Fit))/length(RepPList4Fit);

PolyFit_FPR = polyfit(log10(PredictedFPR(SIndex_m2:SIndex_p2)),...
        log10(RepPList4Fit(SIndex_m2:SIndex_p2)),1);

alphaCutoff = 10^polyval(PolyFit_FPR,log10(pCutoff));

    
    
    