figure;
spRange = [1 4];
xs = abs(medianEffect(:, 2:3));
xs = xs(intint(goodGenes), :);
corrVector = obs_all(goodGenes, 2);
y = corrVector;
labels = {'cis', 'trans'};

cm = cbrewer2('YlGnBu'); 
cm = cm(end:-1:1, :);
thList = 0:.5:max(x);
for j = 1:2
    x = xs(:, j);
    meanStd = nan(2, length(thList));
    stdNonSymmetric = nan(2, length(thList));
    for i = 1:length(thList)-1
        g = find(x > thList(i) & x < thList(i+1));
        meanStd(1, i) = nanmean(corrVector(g));
        meanStd(2, i) = nanstd(corrVector(g));
        stdNonSymmetric(1, i) = nanstd(corrVector(corrVector < meanStd(1, i)));
        stdNonSymmetric(2, i) = nanstd(corrVector(corrVector > meanStd(1, i)));
    end
    subplot(spRange(1), spRange(2), j+1); hold on;
    dscatter(x, y); 
    colormap(cm);
%     shadedErrorBar(thList, meanStd(1, :), [stdNonSymmetric(2, :); stdNonSymmetric(1, :)], '.-b');
    plot(thList, meanStd(1, :), '.-b', 'linewidth', lw, 'markerSize', 10);
%     scatter(thList, meanStd(1, :), [], 'filled', 'b');
    alpha(.5);
    xlabel(['|change in expression levels, ', labels{j}, '|']);
    ylabel('Regulatory similarity, species');
    axis square; box on; ylim([-1 1.2]); grid on;
    set(gca, 'fontsize', fs);
    [r, p] = corr(x, y, 'rows', 'pairwise');
    disp([labels{j}, ': R = ', num2str(round(r, Nround)), ', p-value = ', num2str(round(p, Nround))]);
end

% D
deType = {'sp', 'cis', 'trans'};
Ngenes = 100;
corrth = 0.5;
fractionLowCorr = nan(nC, 4);
for i = 1:nC
    for j = 1:3
        currDE = comparable_samples_de.(deType{j})(:, i);
        currDE = currDE(intint(goodGenes), :);
        [s, sortdidx] = sort(abs(currDE), 'descend');
        topN = sortdidx(1:Ngenes);
        corrOfTopN = corrVector(topN);
        currFraction = length(find(corrOfTopN < corrth))./Ngenes;
        fractionLowCorr(i, j) = currFraction;
    end
    randomGenes = datasample(1:length(goodGenes), Ngenes);
    corrOfRandom = corrVector(randomGenes);
    currFraction = length(find(corrOfRandom < corrth))./Ngenes;
    fractionLowCorr(i, 4) = currFraction;
end

subplot(spRange(1), spRange(2), 4);
edges = 0:.05:0.6;
h1 = histcounts(fractionLowCorr(:, 2), edges);
h2 = histcounts(fractionLowCorr(:, 3), edges);
h3 = histcounts(fractionLowCorr(:, 4), edges);
% bar(edges(1:end-1), [h1; h2;h3]', 'edgeColor', 'none');
plot(edges(1:end-1), [h1; h2; h3]', 'linewidth', lw);
xlabel(['Fraction of genes with Rc < ', num2str(corrth)]);
legend({'cis', 'trans', 'random'});
ylabel('# comparable samples');
set(gca, 'fontsize', fs);
% p24 = kruskalwallis(fractionLowCorr(:, 2:4), [], 'off');
% p23 = kruskalwallis(fractionLowCorr(:, 2:3), [], 'off');
p24 = anova1(fractionLowCorr(:, 2:4), [], 'off');
p23 = anova1(fractionLowCorr(:, 2:3), [], 'off');
disp(['ANOVA cis,trans,random: pvalue = ', num2str(p24)]);
disp(['ANOVA cis,trans: pvalue = ', num2str(p23)]);
xl = xlim; yl = ylim;
% text(mean(xl), mean(yl), ['p-value = ', num2str(p)], 'fontsize', fs);
set(gcf, 'color', 'w', 'position', [1        -116        1920         962]); 
axis square; axis tight; grid on;
