%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab script to process experiments that quantify noise in gene expression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
load Data_exp1.dat
load Data_exp2.dat
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get columns:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
C1 = Data_exp1(:,1);
Y1 = Data_exp1(:,2);
C2 = Data_exp2(:,1);
Y2 = Data_exp2(:,2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot both curves on the same page
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
figure
subplot(2,2,1)
plot(C1,Y1,'or');
axis([0.5 2.5 0 4]);
xlabel('C')
ylabel('Y')
title('Experiment 1')
subplot(2,2,2)
plot(C2,Y2,'or');
axis([0.5 2.5 0 4]);
xlabel('C')
ylabel('Y')
title('Experiment 2')
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute noise levels for each experiment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
cy_mean1 = mean(C1)*mean(Y1);
num1=mean((C1-Y1).^2);
den1 = 2*cy_mean1;
noise_int1 = num1/den1;
num1 = mean(C1.*Y1)-cy_mean1;
den1 = cy_mean1;
noise_ext1 = num1/den1;
num1=mean(C1.^2 + Y1.^2)-2*cy_mean1;
den1 = 2*cy_mean1;
noise_tot1 = num1/den1;
%
cy_mean2 = mean(C2)*mean(Y2);
num1=mean((C2-Y2).^2);
den1 = 2*cy_mean2;
noise_int2 = num1/den1;
num1 = mean(C2.*Y2)-cy_mean2;
den1 = cy_mean2;
noise_ext2 = num1/den1;
num1=mean(C2.^2 + Y2.^2)-2*cy_mean2;
den1 = 2*cy_mean2;
noise_tot2 = num1/den1;