Script2.txt

%This script compares data generated by the model and real data
generation_time = 85;
load('datavariables2.mat'); %This contains both Struan and Simon's data inputs
%load('fourloopgoodshortgene.mat'); % use this if comparing with mutant strains
load('fourloopgood.mat'); %Use this if comparing with WT
chi2_values = zeros(length(av_time_off_input), length(av_time_on_input), length(c_input), length(elongation_input)); %empty vector to be filled
%number_of_cells = 90;%SH9
%number_of_cells = 117;%ADH4
%number_of_cells = 232; %WT
number_of_cells = 100
nucleus = simon_ADH4_nucleus*number_of_cells; %Need to multiply by the number of cells if using Simon's data
cytoplasm = simon_ADH4_cytoplasm*number_of_cells; %Need to multiply by the number of cells if using Simon's data
figure(21)
            axis = [0:(length(nucleus)-1)];
            bar(axis,nucleus,'b')
            xlabel('Number of nuclear transcripts', 'FontName', 'Cambria Math', 'FontSize', 18)
            ylabel('Time (minutes)','FontName', 'Cambria Math', 'FontSize', 18)
            title('Struan WT data nucleus','FontName', 'Cambria Math', 'FontSize', 18)
            xlim([-1 15])
            ylim([0 300])
figure(23)
axis = [0:(length(cytoplasm)-1)];
bar(axis,cytoplasm,'g')
xlabel('Number of cytoplasmic transcripts', 'FontName', 'Cambria Math', 'FontSize', 18)
ylabel('Time (minutes)','FontName', 'Cambria Math', 'FontSize', 18)
title('Struan WT data','FontName', 'Cambria Math', 'FontSize', 18)
xlim([-1,30])
ylim([0 150])
iteration_av_time_off = 1;
for off = 1:length(av_time_off_input);
    iteration_av_time_on = 1;
    for on = 1:length(av_time_on_input);
    iteration_av_time_on ;
    iteration_c_values = 1;
        for c = 1:length(c_input);
            iteration_c_values;
            iteration_retention = 1;
            for elongation = 1:length(elongation_input);
            
%             A = importdata('WT_1.txt');
%             A = A(1:end,1:2);
%             nucleus_real_data = A(1:end,2);
%             cytoplasm_real_data = A(1:end,1);

%THIS IS FOR CONVERTING DATA IF GIVEN IN NON-HISTOGRAM FORMAT. NO NEED IF
%USING DATA I HAVE ALREADY FORMATTED AND SAVED IN datavariables.mat
            
%             F = hist(cytoplasm_real_data, [0:max(cytoplasm_real_data)]);
%             Q = hist(nucleus_real_data, [0:max(nucleus_real_data)]);
%             nucleus_real_data = struan_SH9_nucleus;%*number_of_cells
%             cytoplasm_real_data = struan_SH9_cytoplasm;%*number_of_cells
            
%             F = hist(cytoplasm_real_data, [0:max(cytoplasm_real_data)]);
%             Q = hist(nucleus_real_data, [0:max(nucleus_real_data)]);
%             nucleus = nucleus_real_data;
%             cytoplasm = cytoplasm_real_data;
            
      
            %F = (10/(sum(F)))*hist(cytoplasm_real_data, [0:max(cytoplasm_real_data)]);%have normalised & added *10
            %Q = (10/(sum(Q)))*hist(nucleus_real_data, [0:max(nucleus_real_data)]);
   

            first = iteration_av_time_off
            second = iteration_av_time_on
            third = iteration_c_values
            fourth = iteration_retention
           
            average_nucleus_adjusted = (total_cell_T_array_nucleus{first,second,third,fourth}*((sum(nucleus))/(generation_time)))
            l = abs(length(nucleus)-length(average_nucleus_adjusted));
            if length(average_nucleus_adjusted) < length(nucleus)
                average_nucleus_adjusted = [average_nucleus_adjusted, zeros(1,l)];
            elseif length(average_nucleus_adjusted) > length(nucleus)
                    %Q = [Q, zeros(1,l)]; %Doesn't truncate model
                    average_nucleus_adjusted = average_nucleus_adjusted(1:length(nucleus)); %truncates model
            else 
            end

            average_nucleus_adjusted;
            nucleus;
            % USING A VARIATON OF THE CHI SQUARED TEST, COMPARE THE
            % DISTRIBUTIONS WITH ONE ANOTHER
            number = length(average_nucleus_adjusted);
            chi2 = zeros(1, number);
            for num = 1:number
                if nucleus(num) <=1;
                    chi2(num) = 0; %!!!!!!!!
                else
                chi2(num) = ((average_nucleus_adjusted(num)-nucleus(num))^2)/nucleus(num);
                end
            end
            chi2;
            chi2_sum = sum(chi2); %sums all the transcript number chi2values into one number

            
            v = (length(average_nucleus_adjusted)-1);
            x = chi2_sum;
            q = chi2cdf(x,v);
            p = 1-q;
            chi2_values(iteration_av_time_off, iteration_av_time_on, iteration_c_values, iteration_retention) = p;
            iteration_retention = iteration_retention +1
            end
            iteration_c_values = iteration_c_values+1;
            end
    iteration_av_time_on = iteration_av_time_on+1;
    end
    iteration_av_time_off = iteration_av_time_off+1;
end
chi2_values;
[value, location]  = max(chi2_values(:));
[R,C,S,T] = ind2sub(size(chi2_values),location);
%These are the iterations where you find the best matching distribution
R; %Row (average time off)
C; %Column (average time on)
S; % c iteraton
T; % elongation
figure(22)
            axis = [0:(length(nucleus)-1)];
            bar(axis,total_cell_T_array_nucleus{R,C,S,T}(1:length(nucleus))*((sum(nucleus))/(generation_time)),'b')
          
            xlabel('Number of nuclear transcripts', 'FontName', 'Cambria Math', 'FontSize', 18)
            ylabel('Time (minutes)','FontName', 'Cambria Math', 'FontSize', 18)
            title('Model Wild Type data','FontName', 'Cambria Math', 'FontSize', 18)
            xlim([-1,15])
            ylim([0,300])
            
%REPEAT FOR THE CYTOPLASM

iteration_av_time_off = 1;
for off = av_time_off_input;
    iteration_av_time_off;
    iteration_av_time_on = 1;
    for on = av_time_on_input;
    iteration_av_time_on ;
    iteration_c_values = 1;
        for c = c_input;
            iteration_c_values;
            iteration_retention = 1;
            for elongation = elongation_input;
            
            first = iteration_av_time_off
            second = iteration_av_time_on
            third = iteration_c_values
            fourth = iteration_retention
           
            average_cytoplasm_adjusted = (total_cell_T_array_cytoplasm{first, second,third, fourth}*((sum(cytoplasm))/(generation_time)));
            
l = abs((length(cytoplasm))-length(average_cytoplasm_adjusted));
            if length(average_cytoplasm_adjusted) < length(cytoplasm)
                average_cytoplasm_adjusted = [average_cytoplasm_adjusted, zeros(1,l)];
            elseif length(average_cytoplasm_adjusted) > length(cytoplasm)
                 average_cytoplasm_adjusted = average_cytoplasm_adjusted(1:length(cytoplasm));
            else 
            end

            number = length(average_cytoplasm_adjusted);
            chi2_cytoplasm = zeros(1, number);
            for num = 1:number
                if cytoplasm(num) <=1;
                    chi2_cytoplasm(num) = 0; %!!!!!!!!
                else
                chi2_cytoplasm(num) = ((average_cytoplasm_adjusted(num)-cytoplasm(num))^2)/cytoplasm(num);
                end
            end
chi2_cytoplasm;
            chi2_sum_cytoplasm = sum(chi2_cytoplasm);

            v = (length(average_cytoplasm_adjusted)-1);
            x = chi2_sum_cytoplasm;
            q = chi2cdf(x,v);
            p = 1-q;
            chi2_values_cytoplasm(iteration_av_time_off, iteration_av_time_on, iteration_c_values, iteration_retention) = p;
            iteration_retention = iteration_retention +1;
            end
            iteration_c_values = iteration_c_values+1;
            end
    iteration_av_time_on = iteration_av_time_on+1;
    end
    iteration_av_time_off = iteration_av_time_off+1;
end
chi2_values_cytoplasm;
[value_cytoplasm, location_cytoplasm]  = max(chi2_values_cytoplasm(:));
[R_cytoplasm,C_cytoplasm,S_cytoplasm, T_cytoplasm] = ind2sub(size(chi2_values_cytoplasm),location_cytoplasm);
R_cytoplasm;
C_cytoplasm;
S_cytoplasm;
T_cytoplasm;
figure(24)
            axis = [0:29];
            bar(axis,total_cell_T_array_cytoplasm{R_cytoplasm,C_cytoplasm,S_cytoplasm, T_cytoplasm}(1:30)*((sum(cytoplasm))/(generation_time)),'g')
    xlabel('Number of cytoplasmic transcripts', 'FontName', 'Cambria Math', 'FontSize', 18)
    ylabel('Time (minutes)','FontName', 'Cambria Math', 'FontSize', 18)
    title('Model Wild Type data','FontName', 'Cambria Math', 'FontSize', 18)
    xlim([-1,30])
    ylim([0 150])