% Function: This is the program that will perform the full statistical
% analysis for a set of experiments

% Inputs:
% FileCell - A cell array with the filenames
% WellInfoCell - A nested cell array with the well locations for the 
%    different experimental conditions
% ShiftFormat - A cell array made of lists that rearranges the data from 
%    the xml order to a desired format for display figures 
% NumWellCell - A cell array made of lists that contain the number of wells
%    on each plate
% BeadRegionCell - A cell array made of lists that contain the bead region
%    numbers for all of the IPs
% IPNames - A cell where each entry is a string naming the respective IP
%    target
% ProbeNames - A cell where each entry is a string naming the respective
%    Probe target
% ExpId - A cell where each entry is a string that names each experiment
% CondId - A cell where each entry is a string that names the different
%    experimental conditions tested
% CellType - A string with a general name for the experiment type

% Outputs:
% This program will output the data as matlab figures, csv files, and will
% add variables to the workspace

function [] = MultiplexFullAnalysis_2or3Plates_wXML(FileCell,WellInfoCell,...
    ShiftFormat,NumWellCell,BeadRegionCell,IP_Names,Probe_Names,...
    ExpId,CondId,CellType)

ExpNum = length(FileCell);

% Default Values (For visualizations)
MedianMin = 1;
MedianMax = 20000;
PMin = 0;
PMax = 1;
FoldMin = -2;
FoldMax = 2;

% Values used in each loop
ConditionNum = length(CondId);
CondComb = combnk(1:ConditionNum,2);
[CombinationNum,~] = size(CondComb);

Log2FCCell = cell(1,ExpNum);
PValueCell = cell(1,ExpNum);
RepLog2FCCell = cell(1,ExpNum);
RepPValueCell = cell(1,ExpNum);
MedAvgDataCell = cell(1,ExpNum);

for i = 1:ExpNum
    
    RepLog2FCCelli = cell(ConditionNum,1);
    RepPValueCelli = cell(ConditionNum,1);

    ShiftFormati = ShiftFormat{i};
    
    if length(NumWellCell{i}) == 3
         
        % Define Values
        BeadRegion = BeadRegionCell{i};
        Identifier = ExpId{i};

        % Import Plate Data
        [Plate1_Data] = MultiPlexAnlys_XMLImport(FileCell{i}{1},...
            [5000,25000],BeadRegion); %Function - ImportFormat
        [Plate2_Data] = MultiPlexAnlys_XMLImport(FileCell{i}{2},...
            [5000,25000],BeadRegion); %Function - ImportFormat
        [Plate3_Data] = MultiPlexAnlys_XMLImport(FileCell{i}{3},...
            [5000,25000],BeadRegion); %Function - ImportFormat

        ConcatDatai = cell(ConditionNum,1);
        MedAvgDatai = cell(ConditionNum,1);

        for j = 1:ConditionNum

            [Condj_Med1,Condj_Dist1] = MultiPlexMatrixCold_3Plates(...
                Plate1_Data,Plate2_Data,Plate3_Data,WellInfoCell{i}{j}{1}{1},...
                WellInfoCell{i}{j}{1}{2},WellInfoCell{i}{j}{1}{3},...
                BeadRegion,ShiftFormati{1},ShiftFormati{2},ShiftFormati{3});
            %Function - ImportFormat
            
            [Condj_Med2,Condj_Dist2] = MultiPlexMatrixCold_3Plates(...
                Plate1_Data,Plate2_Data,Plate3_Data,WellInfoCell{i}{j}{2}{1},...
                WellInfoCell{i}{j}{2}{2},WellInfoCell{i}{j}{2}{3},...
                BeadRegion,ShiftFormati{1},ShiftFormati{2},ShiftFormati{3});
            %Function - ImportFormat

            assignin('base', strcat(CondId{j},'_Med1',Identifier),...
                Condj_Med1)
            assignin('base', strcat(CondId{j},'_Dist1',Identifier),...
                Condj_Dist1)
            assignin('base', strcat(CondId{j},'_Med2',Identifier),...
                Condj_Med2)
            assignin('base', strcat(CondId{j},'_Dist2',Identifier),...
                Condj_Dist2)

            Condj_MedAvg = MatrixMean(Condj_Med1,Condj_Med2); 
            %Function - Statistics             
            Condj_Concat = ConcatRep(Condj_Dist1,Condj_Dist2); 
            %Function - ImportFormat        

            assignin('base', strcat(CondId{j},'_MedAvg',Identifier),...
                Condj_MedAvg);
            assignin('base', strcat(CondId{j},'_Concat',Identifier),...
                Condj_Concat);

            MedAvgDatai{j} = Condj_MedAvg;
            ConcatDatai{j} = Condj_Concat;

            [FC_RepCondj] = Fold_Matrix(Condj_Med1,Condj_Med2); 
            %Function - Statistics
            [Log2FC_RepCondj] = log2(FC_RepCondj);
            assignin('base', strcat('Log2FC_', CondId{j},Identifier),...
                Log2FC_RepCondj);

            [PValue_RepCondj] = KS_P_Matrix(Condj_Dist1,Condj_Dist2); 
            %Function - Statistics
            assignin('base', strcat('PValue_Rep', CondId{j},Identifier),...
                PValue_RepCondj);

            RepLog2FCCelli{j} = Log2FC_RepCondj;
            RepPValueCelli{j} = PValue_RepCondj;

            figure(01)
            PlotMatrixINV(Condj_MedAvg,strcat(...
                'Multiplex Median Fluorescence of',{' '},CondId{j},{' '},...
                'for the Experiment with',{' '},CellType,{' '},'on',{' '},...
                Identifier),Probe_Names,IP_Names,...
                0,[MedianMin, MedianMax],'Colormap.txt') 
            %Function - Graphical Displays
            saveas(01,strcat(CondId{j},'_Median',Identifier,'.fig'),'fig'); 
            close(gcf)

%           Commented out code. This code can be used to evaluate the 
%           p-values when comparing replicate samples
%             figure(02)
%             PlotMatrixINV(PValue_RepCondj,strcat(...
%                 'Multiplex Replicate KS P-Values',{' '},CondId{j},{' '},...
%                 'for the Experiment with',{' '},CellType,{' '},'on',{' '},...
%                 Identifier),Probe_Names,IP_Names,1,[PMin, PMax],...
%                 'PMap_441.txt')
%             %Function - Graphical Displays
%             saveas(02,strcat(CondId{j},'RepP',Identifier,'.fig'),'fig');
%             close(gcf)

        end
    elseif length(NumWellCell{i}) == 2 
        
        % Define Values
        BeadRegion = BeadRegionCell{i};
        Identifier = ExpId{i};

        % Import Plate Data
        [Plate1_Data] = MultiPlexAnlys_XMLImport(FileCell{i}{1},...
            [5000,25000],BeadRegion); %Function - ImportFormat 
        [Plate2_Data] = MultiPlexAnlys_XMLImport(FileCell{i}{2},...
            [5000,25000],BeadRegion); %Function - ImportFormat 

        ConcatDatai = cell(ConditionNum,1);
        MedAvgDatai = cell(ConditionNum,1);
        
        for j = 1:ConditionNum

            [Condj_Med1,Condj_Dist1] = MultiPlexMatrixCold(Plate1_Data,...
            Plate2_Data,WellInfoCell{i}{j}{1}{1},...
            WellInfoCell{i}{j}{1}{2},BeadRegion,...
            ShiftFormati{1},ShiftFormati{2}); 
            %Function - ImportFormat

            [Condj_Med2,Condj_Dist2] = MultiPlexMatrixCold(Plate1_Data,...
            Plate2_Data,WellInfoCell{i}{j}{2}{1},...
            WellInfoCell{i}{j}{2}{2},BeadRegion,...
            ShiftFormati{1},ShiftFormati{2});
            %Function - ImportFormat

            assignin('base', strcat(CondId{j},'_Med1',Identifier),...
                Condj_Med1)
            assignin('base', strcat(CondId{j},'_Dist1',Identifier),...
                Condj_Dist1)
            assignin('base', strcat(CondId{j},'_Med2',Identifier),...
                Condj_Med2)
            assignin('base', strcat(CondId{j},'_Dist2',Identifier),...
                Condj_Dist2)

            Condj_MedAvg = MatrixMean(Condj_Med1,Condj_Med2); 
            %Function - Statistics
            Condj_Concat = ConcatRep(Condj_Dist1,Condj_Dist2); 
            %Function - ImportFormat

            assignin('base', strcat(CondId{j},'_MedAvg',Identifier),...
                Condj_MedAvg);
            assignin('base', strcat(CondId{j},'_Concat',Identifier),...
                Condj_Concat);

            MedAvgDatai{j} = Condj_MedAvg;
            ConcatDatai{j} = Condj_Concat;

            [FC_RepCondj] = Fold_Matrix(Condj_Med1,Condj_Med2); 
            %Function - Statistics
            [Log2FC_RepCondj] = log2(FC_RepCondj);
            assignin('base', strcat('Log2FC_', CondId{j},Identifier),...
                Log2FC_RepCondj);

            [PValue_RepCondj] = KS_P_Matrix(Condj_Dist1,Condj_Dist2);
            %Function - Statistics
            assignin('base', strcat('PValue_Rep', CondId{j},Identifier),...
                PValue_RepCondj);

            RepLog2FCCelli{j} = Log2FC_RepCondj;
            RepPValueCelli{j} = PValue_RepCondj;

            figure(01)
            PlotMatrixINV(Condj_MedAvg,strcat(...
                'Multiplex Median Fluorescence of',{' '},CondId{j},{' '},...
                'for the Experiment with',{' '},CellType,{' '},'on',{' '},...
                Identifier),Probe_Names,IP_Names,...
                1,[0,20000],'Colormap.txt')
            %Function - Graphical Displays
            saveas(01,strcat(CondId{j},'_Median',Identifier,'.fig'),'fig'); 
            close(gcf)

%           Commented out code. This code can be used to evaluate the 
%           p-values when comparing replicate samples
%             figure(02)
%             PlotMatrixINV(PValue_RepCondj,strcat(...
%                 'Multiplex Replicate KS P-Values',{' '},CondId{j},{' '},...
%                 'for the Experiment with',{' '},CellType,{' '},'on',{' '},...
%                 Identifier),Probe_Names,IP_Names,1,[PMin, PMax],'PMap_441.txt') 			%GrahicalDisplay
%             %Function - Graphical Displays
%             saveas(02,strcat(CondId{j},'RepP',Identifier,'.fig'),'fig');
%             close(gcf)

        end
    end
    
    Log2FCCelli = cell(CombinationNum,1);
    PValueCelli = cell(CombinationNum,1);
    
    for alpha = 1:CombinationNum

        % PValues and Fold Changes
       
        % PValue Between Conditions
        [PValue_Comba] = KS_P_Matrix(ConcatDatai{CondComb(alpha,2)},...
            ConcatDatai{CondComb(alpha,1)}); 
        %Function - Statistics
  
        assignin('base', strcat('PValue_',CondId{CondComb(alpha,2)},...
            'v',CondId{CondComb(alpha,1)},Identifier),PValue_Comba);
   
        % Fold Change Between Conditions
        [FC_Comba] = Fold_Matrix(MedAvgDatai{CondComb(alpha,2)},...
            MedAvgDatai{CondComb(alpha,1)});
        %Function - Statistics

        assignin('base', strcat('FC_',CondId{CondComb(alpha,2)},...
            'v',CondId{CondComb(alpha,1)},Identifier),FC_Comba);

        % Log2FoldChange Between Conditions
        [Log2FC_Comba] = log2(FC_Comba);

        assignin('base', strcat('Log2FC_',CondId{CondComb(alpha,2)},...
            'v',CondId{CondComb(alpha,1)},Identifier),Log2FC_Comba)
        
        Log2FCCelli{alpha} = Log2FC_Comba;
        PValueCelli{alpha} = PValue_Comba;
  
        % Figures
        figure(03)
        PlotMatrixINV(Log2FC_Comba,strcat(...
            'Multiplex Log 2 Fold Change Between',{' '},...
            CondId{CondComb(alpha,2)},{' '},'and',{' '},...
            CondId{CondComb(alpha,1)},{' '},'for the Experiment with',...
            {' '},CellType,{' '},'on',{' '},Identifier),...
            Probe_Names,IP_Names,1,[FoldMin, FoldMax],'FoldChange.txt')
        %Function - Graphical Displays
        saveas(03,strcat(CondId{CondComb(alpha,2)},'v',...
            CondId{CondComb(alpha,1)},'FC',Identifier,'.fig'),'fig');
        close(gcf)
        
        figure(04)
        PlotMatrixINV(PValue_Comba,strcat(...
            'Multiplex KS P-Values Between',{' '},...
            CondId{CondComb(alpha,2)},{' '},'and',{' '},...
            CondId{CondComb(alpha,1)},{' '},'for the Experiment with',...
            {' '},CellType,{' '},'on',{' '},Identifier),...
            Probe_Names,IP_Names,1,[PMin, PMax],'PMap_441.txt')
        %Function - Graphical Displays
        saveas(04,strcat(CondId{CondComb(alpha,2)},'v',...
            CondId{CondComb(alpha,1)},'P',Identifier,'.fig'),'fig');
        close(gcf)
        
    end
    
    Log2FCCell{i} = Log2FCCelli;
    PValueCell{i} = PValueCelli;
    RepLog2FCCell{i} = RepLog2FCCelli;
    RepPValueCell{i} = RepPValueCelli;
    MedAvgDataCell{i} = MedAvgDatai;
    
end

assignin('base', strcat('Log2FCCell',CellType),Log2FCCell)
assignin('base', strcat('PValueCell',CellType),PValueCell)
assignin('base', strcat('RepLog2FCCell',CellType),RepLog2FCCell)
assignin('base', strcat('RepPValueCell',CellType),RepPValueCell)
assignin('base', strcat('MedAvgDataCell',CellType),MedAvgDataCell)

IPProbeNum = length(IP_Names)*length(Probe_Names);
NRepValues = ExpNum*ConditionNum*IPProbeNum;

[EmpericalAlphaCell,PerAgreeCell,RepList4Fit,~] = ...
    EmpAdjPValueList(0.05/IPProbeNum,RepPValueCell,NRepValues); 
%Function - Statistics/MetaConsistency

assignin('base', strcat('EmpericalAlphaCell',CellType),EmpericalAlphaCell)
assignin('base', strcat('PerAgreeCell',CellType),PerAgreeCell)
assignin('base', strcat('RepList4Fit',CellType),RepList4Fit)

% Commented out code. This code can be used to troubleshoot and identify 
% probes with consistent variation between replicates.
% assignin('base', strcat('BadProbeCell',CellType),BadProbeCell)

CombId = cell(CombinationNum,1);
for beta = 1:CombinationNum
    CombId{beta} = strcat(CondId{CondComb(beta,2)},'v',...
            CondId{CondComb(beta,1)});
end

IPProbeList = LabelListFormat(Probe_Names,IP_Names);
%Function - ImportFormat

[ExpHitCells] = GlobalConsistentHits(PValueCell,Log2FCCell,...
    EmpericalAlphaCell,PerAgreeCell,IPProbeList,CombId,CellType); 
%Function - Statistics/MetaConsistency
assignin('base', strcat('ExpHitCells',CellType),ExpHitCells)

% Commented out code. This code can be used to troubleshoot and identify 
% probes with consistent variation between replicates.
% [RepHitCells] = GlobalConsistentHits(RepPValueCell,RepLog2FCCell,...
%     EmpericalAlphaCell,PerAgreeCell,IPProbeList,CondId,...
%     strcat(CellType,'Reps'));
% %Function - Statistics/MetaConsistency
% assignin('base', strcat('RepHitCells',CellType),RepHitCells)

% Define RepHitCells as an empty list as it is called in future functions.
RepHitCells = {};

% Generate CSV output
ExcelHitOutput(ExpHitCells,RepHitCells,PValueCell,Log2FCCell,...
    ExpNum,CombinationNum,IPProbeNum,CombId,strcat(CellType,'_Hits.csv'));
% Function - ExcelOutput

ExcelAllOutput(ExpHitCells,PValueCell,Log2FCCell,...
    ExpNum,CombinationNum,IPProbeNum,CombId,IPProbeList,...
    strcat(CellType,'All.csv'));
% Function - ExcelOutput

ExcelMFIOutput(MedAvgDataCell,ExpNum,ConditionNum,IPProbeNum,...
    CondId,IPProbeList,strcat(CellType,'_MFI.csv'));
% Function - ExcelOutput

AvgLog2FCCell = cell(CombinationNum,1);
[ProbeNum,IPNum] = size(Log2FCCell{1}{1}); 

% Generate Avergae FC Matrix
for gamma = 1:CombinationNum
        
        Log2FCAvg_gamma = zeros(ProbeNum,IPNum); 
        
        Log2FCList_gamma = zeros(ExpNum,ProbeNum*IPNum);
        
        for delta = 1:ExpNum
            Log2FCList_gamma(delta,:) = Log2FCCell{delta}{gamma}(:);
        end
        
        Log2FCAvg_gamma(:) = nanmean(Log2FCList_gamma);
        
        AvgLog2FCCell{gamma} = Log2FCAvg_gamma;
        
end
assignin('base', strcat('AvgLog2FCCell',CellType),AvgLog2FCCell)