% figure 2
%% B,C

geneA = 'RPL1A';
geneB = 'RPL1B';
gnA = GP.gene_table.(geneA);
gnB = GP.gene_table.(geneB);

spLocs = {1, 4, [2,3,5,6]};
spRange = [2,3];
fs = 10;
cm1 = cbrewer2('Greys'); cm1 = cm1(end:-1:1, :);
cm2 = cbrewer2('YlGnBu'); cm2 = cm2(end:-1:1, :);
figure;
for i = 1:2
    d1 = expStruct.(gb{i}).exp(gnA, :)';
    d2 = expStruct.(gb{i}).exp(gnB, :)';
    subplot(spRange(1), spRange(2), spLocs{i});
    dscatter(d1, d2);
    colormap(cm1);
    xlabel(geneA);
    ylabel(geneB);
    r = corr(d1, d2, 'rows', 'pairwise');
    title([gb{i}, ', R = ', num2str(round(r, 2))]);
    axis square; grid on; box on;
    set(gca, 'fontsize', fs);
    cb = colorbar; cb.Label.String = 'Density of samples'; cb.Ticks = []; cb.FontSize = fs;
end
   
gnA1 = find(intint==gnA);
gnB1 = find(intint==gnB);
d1 = PCMs.WCDcer.PCM(:, gnA1);
d2 = PCMs.WCDpar.PCM(:, gnA1);
s3 = subplot(spRange(1), spRange(2), spLocs{3});
dscatter(d1, d2); hold on;
scatter(d1(gnB1), d2(gnB1), [], 'r', 'filled');
text(d1(gnB1), d2(gnB1), geneB, 'fontsize', fs);
colormap(s3, cm2);
set(gca,'xtick', -1:.5:1, 'ytick', -1:.5:1);
axis square; grid on; box on;
axis([-1 1 -1 1]);
xlabel('correlation, cer'); ylabel('correlation, par');
r = corr(d1, d2, 'rows', 'pairwise');
title(['Regulatory similarity of ', geneA, ', R = ', num2str(round(r, 2))]);
set(gca, 'fontsize', fs);
set(gcf, 'color', 'w', 'position', [520   369   891   429]);
cb = colorbar; cb.Label.String = 'Density of genes'; cb.Ticks = [];cb.FontSize = fs;
%% D

toPrint = 'parents';
rand1 = datasample(1:length(goodGenes), length(goodGenes), 'Replace', false);
switch toPrint
    case 'parents'
        cerRand = PCMs.WCDcer.PCM(rand1, :);
        c = corr(cerRand, PCMs.WCDpar.PCM(goodGenes, :), 'rows', 'pairwise');
        drandomPairs = diag(c);
        nn = nanmedian(nnWithinGenome.secondOrder.corrs(:, 1:2), 2);
%         dists = [obs_all(:, 1), nn, expP, drandomPairs];
        dists = [obs_all(:, 1), expP, nn ];
        dists = dists(goodGenes, :);
        legend1 = 'between species';
    case 'hybrid'
        hycRand = PCMs.WCDhyc.PCM(rand1, :);
        c = corr(hycRand, PCMs.WCDhyp.PCM(goodGenes, :), 'rows', 'pairwise');
        drandomPairs = diag(c);
        nn = nanmedian(nnWithinGenome.secondOrder.corrs(:, 3:4), 2);
%         dists = [obs_all(:, 2), nn, expH, drandomPairs];
        dists = [obs_all(:, 2),expH, nn];
        dists = dists(goodGenes, :);
        legend1 = 'between alleles';
end

figure; hold on;
for i = 1:size(dists, 2)
    d = dists(:, i);
    [ycdf, xcdf] = cdfcalc(d);
    yccdf = ycdf(1:end-1);
    plot(xcdf, yccdf, '.-');
end
xlim([-.8 1]);
plot(xlim, [.5 .5], 'k');
set(gca, 'ytick', 0:0.25:1);
legend(horzcat(legend1, {'dataset control', 'nearest neighbor'}), 'location', 'northwest');
xlabel('correlation (Rc)'); ylabel('cumulative fraction of genes');
grid on; axis square; box on; set(gcf, 'color', 'w');  
set(gca, 'fontsize', fs);

th = 0.5;
sum(dists < th, 1)./size(dists,1)

%% E
x = expP;
y = obs_all(:, 1);
figure; dscatter(x, y); hold on;
colormap(cm2);
plot([.2 1], [.2 1], 'k');
set(gca, 'xtick', .2:.2:1, 'ytick', -1:.5:1);
axis square; box on; grid on;
set(gcf, 'color', 'w');
xlabel('dataset control');
ylabel('between species');
cb = colorbar; cb.Label.String = 'Density of genes'; cb.Ticks = [];cb.FontSize = fs;
set(gca, 'fontsize', fs);

th = 0.2;
idx1 = find(y > th);
idx2 = find(y < th);
x1 = x(idx1);
y1 = y(idx1);
[r, p] = corr(x1, y1, 'rows', 'pairwise')

x2 = x(idx2);
y2 = y(idx2);
[r, p] = corr(x2, y2, 'rows', 'pairwise')

% p1 = linortfit2(x1,y1);
% xs = [.5 1];
% f1 = polyval(p1, xs);
% plot(xs, f1, 'r--')
% 
% p2 = linortfit2(x2,y2);
% xs = [.5 1];
% f2 = polyval(p2, xs);
% plot(xs, f2, 'b--')
% 
% legend({'genes', 'x=y', 'fit for y>0.2', 'fit for y<0.2'});

