% STEP 5: Statistics Tabulation
%
% The focus of Step 5 is to perform the relevant statistics and generate other information for each of the 1023 possible isoforms
%
% The following data will be reported here via an output file.
%
% A. Decode base 10 values (1..1023) to exon signatures.
% B. Counts for each of the 1023 possible isoforms for each of the eight slides.
% C. Z-scores and P-values for comparison of summed isoform counts for polonies derived from transfected vs. untransfected cell lines.
% D. Counts for individual exons for each of the eight slides
% E. Z-scores and P-values for comparison of summed exon counts for polonies derived from transfected vs. untransfected cell lines.

% 1. Reset variables & close image windows to free up memory.

clear, close all;

% 2. Decode base 10 values (1..1023) to exon signatures.  The goal here is to simply create a cell-array of strings that describe the exon-signature of each isoform
%    that will be used later in this program.

SIG_STRINGS = cell(1023,1);
SIG_BITS = zeros(1023,10);

for i = 1:1023

    T1 = i;
    
    T2 = zeros(10,1);
    
    for j = 10:-1:1    
        
        T3 = 2^(j-1);
        
        if (T1 >= T3)
            T2(j) = 1;
            SIG_BITS(i,j) = 1;
            T1 = T1 - T3;
        end        
        
    end
    
    S = '';
    
    for j = 1:10
    
        if (T2(j) == 1)
            
            S = strcat(S,num2str(j));    
            
        else
            
            S = strcat(S,'-');
            
        end
        
        if (j < 10)
        
            S = strcat(S,'-');    
            
        end 
        
    end
    
    SIG_STRINGS{i} = S;

end

% 3. Load list of called polonies

load ('final_polony_list.mat');

POLONY_LIST = FINAL_POLONY_LIST;

% 4. Now go back and fill in how many polonies derived from each of the ten most commonly observed isoforms was observed on each of the eight slides

COUNTS = zeros(1023,8);

for (i = 1:1023)
    for (j = 1:8)

        C = size(find(POLONY_LIST(:,1) == j & POLONY_LIST(:,2) == i));
        COUNTS(i,j) = C(1);
                
    end
end

% 5. Calculate confidence interval of fold change (Z-scores and P-values)

COL_A = COUNTS(:,1) + COUNTS(:,3) + COUNTS(:,5) + COUNTS(:,7); % Sum counts from slides 1 & 3
COL_B = COUNTS(:,2) + COUNTS(:,4) + COUNTS(:,6) + COUNTS(:,8); % Sum counts from slides 2 & 4

M = sum(COL_A);
N = sum(COL_B);

P1 = COL_A / M;
P2 = COL_B / N;
P = ((M*P1)/(M+N)) + ((N*P2)/(M+N));
Q = 1 - P;

Z = (P1 - P2)./sqrt(P.*Q*((1/M) + (1/N)));

PV = normcdf(-1*abs(Z));

% 6. Add this information to the array of counts

I_STATS = horzcat(COUNTS,Z,PV);

% 7. Now start counts of indiviudal exons

E_STATS = zeros(10,8);

for (i = 1:1023)

    for (j = 1:10)
    
        if (SIG_BITS(i,j) == 1)
            
            E_STATS(j,:) = E_STATS(j,:) + I_STATS(i,1:8);
            
        end
        
    end
    
end

% 8. Calculate confidence interval of fold change of exon counts (Z-scores and P-values)

COL_A = E_STATS(:,1) + E_STATS(:,3) + E_STATS(:,5) + E_STATS(:,7); % Sum counts from slides 1 & 3
COL_B = E_STATS(:,2) + E_STATS(:,4) + E_STATS(:,6) + E_STATS(:,8); % Sum counts from slides 2 & 4

M = sum(COL_A);
N = sum(COL_B);

P1 = COL_A / M;
P2 = COL_B / N;
P = ((M*P1)/(M+N)) + ((N*P2)/(M+N));
Q = 1 - P;

Z = (P1 - P2)./sqrt(P.*Q*((1/M) + (1/N)));

PV = normcdf(-1*abs(Z));

% 9. Add this information to the array of exon counts

E_STATS = horzcat(E_STATS,Z,PV);

% 10. Save important variables in both MATLAB and ASCII forms

save('final_stats.mat','I_STATS','E_STATS','SIG_STRINGS','SIG_BITS');
save('isoform_stats.txt','I_STATS','-ASCII','-DOUBLE','-TABS');
save('exon_stats.txt','E_STATS','-ASCII','-DOUBLE','-TABS');
save('isoform_bits.txt','SIG_BITS','-ASCII','-DOUBLE','-TABS');

fid = fopen('isoform_strings.txt','w');
for i = 1:1023
    fprintf(fid,'%s\n',SIG_STRINGS{i});
end
fclose(fid);

% 11. Reset variables & close image windows to free up memory.

clear, close all

 

 