%Import data

RB_Z1_Foxp2 = [RB1Foxp2;RB3Foxp2;RB4Foxp2];
RB_Z2_Foxp2 = [RB5Foxp2;RB6Foxp2;RB8Foxp2];
RB_Z3_Foxp2 = [RB9Foxp2;RB12Foxp2;RB14Foxp2];

RB_Z1_Pou6f2 = [RB1Pou6f2;RB3Pou6f2;RB4Pou6f2];
RB_Z2_Pou6f2 = [RB5Pou6f2;RB6Pou6f2;RB8Pou6f2];
RB_Z3_Pou6f2 = [RB9Pou6f2;RB12Pou6f2;RB14Pou6f2];

RF_Z1_Foxp2 = [RF2Foxp2;RF3Foxp2;RF4Foxp2];
RF_Z2_Foxp2 = [RF5Foxp2;RF7Foxp2;RF8Foxp2];
RF_Z3_Foxp2 = [RF9Foxp2;RF12Foxp2;RF13Foxp2];

RF_Z1_Pou6f2 = [RF2Pou6f2;RF3Pou6f2;RF4Pou6f2];
RF_Z2_Pou6f2 = [RF5Pou6f2;RF7Pou6f2;RF8Pou6f2];
RF_Z3_Pou6f2 = [RF9Pou6f2;RF12Pou6f2;RF13Pou6f2];

A_Z1_Foxp2 = [A5Foxp2;A6Foxp2;A7Foxp2];
A_Z2_Foxp2 = [A8Foxp2;A9Foxp2;A10Foxp2];
A_Z3_Foxp2 = [A11Foxp2;A12Foxp2;A13Foxp2];

A_Z1_Pou6f2 = [A5Pou6f2;A6Pou6f2;A7Pou6f2];
A_Z2_Pou6f2 = [A8Pou6f2;A9Pou6f2;A10Pou6f2];
A_Z3_Pou6f2 = [A11Pou6f2;A12Pou6f2;A13Pou6f2];

% Aggregate data across brains
Z1_Foxp2 = [RB_Z1_Foxp2;RF_Z1_Foxp2;A_Z1_Foxp2];
Z2_Foxp2 = [RB_Z2_Foxp2;RF_Z2_Foxp2;A_Z2_Foxp2];
Z3_Foxp2 = [RB_Z3_Foxp2;RF_Z3_Foxp2;A_Z3_Foxp2];

Z1_Pou6f2 = [RB_Z1_Pou6f2;RF_Z1_Pou6f2;A_Z1_Pou6f2];
Z2_Pou6f2 = [RB_Z2_Pou6f2;RF_Z2_Pou6f2;A_Z2_Pou6f2];
Z3_Pou6f2 = [RB_Z3_Pou6f2;RF_Z3_Pou6f2;A_Z3_Pou6f2];


%

% Plot by brain

%Plot Z1_Foxp2 by brain
figure;
a = scatter(RB_Z1_Foxp2(:,1),RB_Z1_Foxp2(:,2),10,'green','filled');
hold on
b = scatter(RF_Z1_Foxp2(:,1),RF_Z1_Foxp2(:,2),10,'blue','filled');
c = scatter(A_Z1_Foxp2(:,1),A_Z1_Foxp2(:,2),10,'red','filled');
legend('RB','RF','A')
title('Foxp2 Z1 by brain')
xlim([500 2500]);
ylim([-2500 -500]);

%Plot Z2_Foxp2 by brain
figure;
a = scatter(RB_Z2_Foxp2(:,1),RB_Z2_Foxp2(:,2),10,'green','filled');
hold on
b = scatter(RF_Z2_Foxp2(:,1),RF_Z2_Foxp2(:,2),10,'blue','filled');
c = scatter(A_Z2_Foxp2(:,1),A_Z2_Foxp2(:,2),10,'red','filled');
legend('RB','RF','A')
title('Foxp2 Z2 by brain')
xlim([500 2500]);
ylim([-2500 -500]);

%Plot Z3_Foxp2 by brain
figure;
a = scatter(RB_Z3_Foxp2(:,1),RB_Z3_Foxp2(:,2),10,'green','filled');
hold on
b = scatter(RF_Z3_Foxp2(:,1),RF_Z3_Foxp2(:,2),10,'blue','filled');
c = scatter(A_Z3_Foxp2(:,1),A_Z3_Foxp2(:,2),10,'red','filled');
legend('RB','RF','A')
title('Foxp2 Z3 by brain')
xlim([500 2500]);
ylim([-2500 -500]);

%Plot Z1_Pou6f2 by brain
figure;
a = scatter(RB_Z1_Pou6f2(:,1),RB_Z1_Pou6f2(:,2),10,'green','filled');
hold on
b = scatter(RF_Z1_Pou6f2(:,1),RF_Z1_Pou6f2(:,2),10,'blue','filled');
c = scatter(A_Z1_Pou6f2(:,1),A_Z1_Pou6f2(:,2),10,'red','filled');
legend('RB','RF','A')
title('Pou6f2 Z1 by brain')
xlim([500 2500]);
ylim([-2500 -500]);

%Plot Z2_Pou6f2 by brain
figure;
a = scatter(RB_Z2_Pou6f2(:,1),RB_Z2_Pou6f2(:,2),10,'green','filled');
hold on
b = scatter(RF_Z2_Pou6f2(:,1),RF_Z2_Pou6f2(:,2),10,'blue','filled');
c = scatter(A_Z2_Pou6f2(:,1),A_Z2_Pou6f2(:,2),10,'red','filled');
legend('RB','RF','A')
title('Pou6f2 Z2 by brain')
xlim([500 2500]);
ylim([-2500 -500]);

%Plot Z3_Pou6f2 by brain
figure;
a = scatter(RB_Z3_Pou6f2(:,1),RB_Z3_Pou6f2(:,2),10,'green','filled');
hold on
b = scatter(RF_Z3_Pou6f2(:,1),RF_Z3_Pou6f2(:,2),10,'blue','filled');
c = scatter(A_Z3_Pou6f2(:,1),A_Z3_Pou6f2(:,2),10,'red','filled');
legend('RB','RF','A')
title('Pou6f2 Z3 by brain')
xlim([500 2500]);
ylim([-2500 -500]);


%

%Plot Foxp2 vs Pou6f2

newgreen = '#07A815';

%Z1
figure;
a = scatter(Z1_Foxp2(:,1),Z1_Foxp2(:,2),10,'filled');
a.MarkerFaceColor=	newgreen;
hold on
b = scatter(Z1_Pou6f2(:,1),Z1_Pou6f2(:,2),10,'magenta','filled');
legend('Foxp2','Pou6f2')
xlim([500 2500]);
ylim([-2500 -500]);

%Z2
figure;
a = scatter(Z2_Foxp2(:,1),Z2_Foxp2(:,2),10,'filled');
a.MarkerFaceColor=	newgreen;
hold on
b = scatter(Z2_Pou6f2(:,1),Z2_Pou6f2(:,2),10,'magenta','filled');
legend('Foxp2','Pou6f2')
xlim([500 2500]);
ylim([-2500 -500]);

%Z3
figure;
a = scatter(Z3_Foxp2(:,1),Z3_Foxp2(:,2),10,'filled');
a.MarkerFaceColor=	newgreen;
hold on
b = scatter(Z3_Pou6f2(:,1),Z3_Pou6f2(:,2),10,'magenta','filled');
legend('Foxp2','Pou6f2')
xlim([500 2500]);
ylim([-2500 -500]);


%

%Calculate density

N_LEVELS = 4;   % how many contour levels should we draw?
N_BINS = 256;   % how big of an NxN matrix to use for kde?

xBounds = linspace(500,2500,N_BINS);
yBounds = linspace(-2500,-500,N_BINS);

[~, density_Z1_Foxp2]=kde2d(Z1_Foxp2(:,(1:2)),N_BINS,[500 -2500],[2500 -500]);
[~, density_Z1_Pou6f2]=kde2d(Z1_Pou6f2(:,(1:2)),N_BINS,[500 -2500],[2500 -500]);

[~, density_Z2_Foxp2]=kde2d(Z2_Foxp2(:,(1:2)),N_BINS,[500 -2500],[2500 -500]);
[~, density_Z2_Pou6f2]=kde2d(Z2_Pou6f2(:,(1:2)),N_BINS,[500 -2500],[2500 -500]);

[~, density_Z3_Foxp2]=kde2d(Z3_Foxp2(:,(1:2)),N_BINS,[500 -2500],[2500 -500]);
[~, density_Z3_Pou6f2]=kde2d(Z3_Pou6f2(:,(1:2)),N_BINS,[500 -2500],[2500 -500]);

%Calculate density difference

densityDifference_Z1 = density_Z1_Foxp2 - density_Z1_Pou6f2;
densityDifference_Z2 = density_Z2_Foxp2 - density_Z2_Pou6f2;
densityDifference_Z3 = density_Z3_Foxp2 - density_Z3_Pou6f2;

%

%Contour plots
    
%Z1 contour plot
figure; set(gcf,'Color','w');
    cax = subplot(1,1,1);   
[~,f] = contourf(xBounds,yBounds,density_Z1_Foxp2,N_LEVELS);
f.Fill = 'off';
f.EdgeColor = 'flat';
f.LineWidth = 1;
f.LineColor = newgreen;     
 hold on   
[~,p] = contourf(xBounds,yBounds,density_Z1_Pou6f2,N_LEVELS);
p.Fill = 'off';
p.EdgeColor = 'flat';
p.LineWidth = 1;
p.LineColor = 'magenta';


%Z2 contour plot
figure; set(gcf,'Color','w');
    cax = subplot(1,1,1);   
[~,f] = contourf(xBounds,yBounds,density_Z2_Foxp2,N_LEVELS);
f.Fill = 'off';
f.EdgeColor = 'flat';
f.LineWidth = 1;
f.LineColor = newgreen;     
 hold on   
[~,p] = contourf(xBounds,yBounds,density_Z2_Pou6f2,N_LEVELS);
p.Fill = 'off';
p.EdgeColor = 'flat';
p.LineWidth = 1;
p.LineColor = 'magenta';

%Z3 contour plot
figure; set(gcf,'Color','w');
    cax = subplot(1,1,1);   
[~,f] = contourf(xBounds,yBounds,density_Z3_Foxp2,N_LEVELS);
f.Fill = 'off';
f.EdgeColor = 'flat';
f.LineWidth = 1;
f.LineColor = newgreen;     
 hold on   
[~,p] = contourf(xBounds,yBounds,density_Z3_Pou6f2,N_LEVELS);
p.Fill = 'off';
p.EdgeColor = 'flat';
p.LineWidth = 1;
p.LineColor = 'magenta';

%

%Density difference plots

%Density difference Z1
figure; set(gcf,'Color','w');
cax = subplot(1,1,1);
[~,j] = contourf(xBounds,yBounds,densityDifference_Z1, N_LEVELS);
j.Fill = 'on';
j.EdgeColor = 'flat';
j.LineWidth = 1;
hcb=colorbar;
set(hcb,'YTick',[])
xlim([500 2500])
ylim([-2500 -500])
j.Fill = 'off';
j.LineWidth=3;
colormap(othercolor('PiYG11'));
caxis([-3e-06 3e-06])

%Density difference Z2
figure; set(gcf,'Color','w');
cax = subplot(1,1,1);
[~,j] = contourf(xBounds,yBounds,densityDifference_Z2, N_LEVELS);
j.Fill = 'on';
j.EdgeColor = 'flat';
j.LineWidth = 1;
hcb=colorbar;
set(hcb,'YTick',[])
xlim([500 2500])
ylim([-2500 -500])
j.Fill = 'off';
j.LineWidth=3;
colormap(othercolor('PiYG11'));
caxis([-3e-06 3e-06])

%Density difference Z3
figure; set(gcf,'Color','w');
cax = subplot(1,1,1);
[~,j] = contourf(xBounds,yBounds,densityDifference_Z3, N_LEVELS);
j.Fill = 'on';
j.EdgeColor = 'flat';
j.LineWidth = 1;
hcb=colorbar;
%set(hcb,'YTick',[])
xlim([500 2500])
ylim([-2500 -500])
j.Fill = 'off';
j.LineWidth=3;
colormap(othercolor('PiYG11'));
caxis([-2e-06 2e-06])



  

%


%Stats Foxp2 vs Pou6f2

%Stats on combined brain (subsample, corr2)