% figure 7 revision
examples_4_scatters('LSO2', {});
expSamples2('LSO2', 'whichModule', 'rb');
plot_enrichment_allSpecies('LSO2', {'Protein synthesis', 'amino-acid biosynthesis'});
%%
gn = find(intint==GP.gene_table.LSO2);
ribo = find(ismember(intint, OGtable.vals{150}));
for i  =1:4
    currEntry = ['WCD', gb{i}];
    currLso2 = PCMs.(currEntry).PCM(gn, :);
    currRibo = PCMs.(currEntry).PCM(ribo, :);
    currRibo = nanmedian(currRibo);
    r = corr(currLso2', currRibo', 'rows', 'pairwise');
    disp(['Lso2 Rc with ribo in ', gb{i}, ' = ', num2str(round(r, 2))]);
  
    currLso2 = dataS.(currEntry)(gn, :);
    currRibo = dataS.(currEntry)(ribo, :);
    currRibo = nanmedian(currRibo);
    r = corr(currLso2', currRibo', 'rows', 'pairwise');
    disp(['Lso2 Re with ribo in ', gb{i}, ' = ', num2str(round(r, 2))]);
end
%% plot along conditions
% gn = GP.gene_table.LSO2;
% dcer = dataS.WCDhyc(gn, :);
% dpar = dataS.WCDhyp(gn, :);
% rib = OGtable.vals{150};
% aa = OGtable.vals{129};
% drib = nanmedian(dataS.WCDhyc(rib, :), 1);
% drib2 = nanmedian(dataS.WCDhyp(rib, :), 1);
% daa = nanmedian(dataS.WCDhyc(aa, :), 1);
% vectors = [dcer; dpar; drib; daa];
% vectors = vectors - nanmedian(vectors, 2);
% [~, s] = sort(drib);
% mm = 50;
% vectors = vectors(:, s);
% for i  =1:size(vectors,1)
%     vectors(i, :) = smoothdata(vectors(i, :), 'movmean', mm);
% end
% colors = [rgb('Blue'); rgb('Purple'); rgb('Green'); rgb('Red')];
% figure; hold on;
% for i  =1:size(vectors,1)
%     plot(vectors(i, :), 'color', colors(i, :));
% end
% axis tight; box on; axis square; grid on;
% set(gcf, 'color', 'w');
% set(gca, 'ytick', -3:1:1);
% ylabel('log2 fold change');
% legend({'LSO2 hyb.cer', 'LSO2 hyb.par', 'Ribosomes', 'AA biosynthesis'}, 'location', 'southwest');
% set(gca, 'fontsize', fs);
%%
exampleGenes = {'LSO2'};
listOfGenes = allGenes(intint(goodGenes));
exampleGenesIdx = ismember_smart(listOfGenes, exampleGenes);

obs_all_gg = obs_all(goodGenes, :);
exp_all_gg = exp_all(goodGenes, :);
load('Rc_genes.mat', 'Rc_cis');
cisGenes = Rc_genes.Rc_cis;
gList = find(ismember(intint(goodGenes), cisGenes));
sz = 50;
fs = 14;
%%
d = obs_all(:, 1) - obs_all(:, 2);
d = d(goodGenes);
cisGenesRc = d(gList);
allOthers = d;
allOthers(gList) = [];
l = 1.5;
edges = -l:.1:l;
h1 = histcounts(allOthers, edges, 'normalization', 'probability');
h2 = histcounts(cisGenesRc, edges, 'normalization', 'probability');
th = .2;
% h3 = histcounts(d(transGenes), edges, 'normalization', 'probability');
figure; hold on;
yl = [0 .4];
% patch([-l -th -th -l], [0 0 yl(2) yl(2)], rgb('tan'));
% patch([th l l th], [0 0 yl(2) yl(2)], rgb('PaleGreen'));
h = bar(edges(1:end-1), [h1; h2]', 'edgeColor', 'none'); 
% plot(edges(1:end-1), [h1; h2]);
xlim([-l l]);	
legend([h(1), h(2)], {'all genes', 'cis-diverging'});
xlabel('Rc species - Rc hybrid');
ylabel('Fraction of Genes');
set(gca, 'fontsize', fs);
axis square; box on; set(gcf, 'color', 'w'); grid on;

th=0.2;
r1 = length(find(abs(allOthers) > th))./length(d);
r2 = length(find(abs(cisGenesRc) > th))./length(gList);
r3 = length(find(allOthers > th))./length(d);
r4 = length(find(cisGenesRc > th))./length(gList);
disp(['all genes > |', num2str(th), '| : ', num2str(r1)]);
disp(['cis genes > |', num2str(th), '| : ', num2str(r2)]);
disp(['all genes > ', num2str(th), ' : ', num2str(r3)]);
disp(['cis genes > ', num2str(th), ' : ', num2str(r4)]);

%% C
x = obs_all_gg(:, 3); % cer-hyc
y = obs_all_gg(:, 4); % par-hyp
z = obs_all_gg(:, 2); % hybrid
% x = zscore(x); y = zscore(y); z = zscore(z);

% there are more cer-like than par-like effects
diffxy = x(gList) - y(gList);
thList = 0:0.1:0.5;
% thList = 0:.5:3;
for i = 1:length(thList)
    currTh = thList(i);
    n1 = length(find(x(gList) - y(gList) > currTh));
    n2 = length(find(y(gList) - x(gList) > currTh));
    disp(['th = ', num2str(currTh), ', x > y : ', num2str(n1), ', x < y : ', num2str(n2)]);
end

figure; 
subplot(1,3,1);
scatter(x, y, sz, [.7 .7 .7]);
hold on;
scatter(x(gList), y(gList), sz, z(gList), 'filled', 'MarkerEdgeColor', 'k');
diffTh = 0.2;
% xl = [min([min(x) min(y)]) 1];
xl = [-0.8 1];
plot(xl, xl + diffTh, 'k--');
plot(xl+diffTh, xl, 'k--');
plot(xl, xl, 'k');
axis([xl(1) xl(2) xl(1) xl(2)]);
cm = cbrewer2('Blues'); cm = cm(end:-1:1, :);
colormap(cm);

xlabel('Rc, cer - hybrid cer');
ylabel('Rc, par - hybrid par');
axis square; grid on; box on; set(gca, 'fontsize', fs);
set(gca, 'xtick', -.5:.5:1, 'ytick', -.5:.5:1);
text(x(exampleGenesIdx), y(exampleGenesIdx), listOfGenes(exampleGenesIdx), 'fontsize', fs);
% text(x(gList), y(gList), listOfGenes(gList));
%% D

nnCorrCer = nnWithinGenome.secondOrder.corrs(:, 3);
nnCorrPar = nnWithinGenome.secondOrder.corrs(:, 4);
x = nnCorrCer(goodGenes);
y = nnCorrPar(goodGenes);

subplot(1,3,2);
scatter(x, y, sz, [.7 .7 .7]);
hold on;
scatter(x(gList), y(gList), sz, z(gList), 'filled', 'MarkerEdgeColor', 'k');
xlabel('Rc, hyb.cer - NNhyb.cer');
ylabel('Rc, hyb.par - NNhyb.par');
axis square; grid on; box on; set(gca, 'fontsize', fs);
text(x(exampleGenesIdx), y(exampleGenesIdx), listOfGenes(exampleGenesIdx), 'fontsize', fs);
axis([.5 1 .5 1]);
set(gca, 'xtick', .5:.1:1, 'ytick', .5:.1:1);
nnTh = 0.85;
plot([nnTh nnTh], ylim, 'k--');
plot(xlim, [nnTh nnTh], 'k--');
% text(x(gList), y(gList), listOfGenes(gList));
%% NN > th
nnTh = 0.85;
x = nnCorrCer(goodGenes);
y = nnCorrPar(goodGenes);
xgty = find(x(gList) > nnTh & y(gList) < nnTh);
ygtx = find(y(gList) > nnTh & x(gList) < nnTh);
n = union(find(x < nnTh), find(y < nnTh));
n2 = union(find(x(gList) < nnTh), find(y(gList) < nnTh));
disp(['# genes < ', num2str(nnTh), ' : ', num2str(length(n)), ', ', num2str(length(n)*100/length(goodGenes)), '%']);
disp(['# genes in gList < ', num2str(nnTh), ' : ', num2str(length(n2)), ', ', num2str(length(n2)*100/length(gList)), '%']);
disp(['# genes x > th, y < th :', num2str(length(xgty))]);
disp(['# genes y > th, x < th :', num2str(length(ygtx))]);

x = obs_all_gg(:, 3); % cer-hyc
y = obs_all_gg(:, 4); % par-hyp
cerLike1 = find(x(gList) - y(gList) > diffTh);
parLike1 = find(y(gList) - x(gList) > diffTh);
nCerlike = intersect(xgty, cerLike1);
nParlike = intersect(ygtx, parLike1);
disp(['# genes x > th, y < th & cer-like: ', num2str(length(nCerlike)), ', ', num2str(length(nCerlike)*100/length(xgty)), '%']);
disp(['# genes y > th, x < th & par-like: ', num2str(length(nParlike)), ', ', num2str(length(nParlike)*100/length(ygtx)), '%']);
%% E

cerByParN = NNtop5.corr(:, 3);
parByCerN = NNtop5.corr(:, 4);
x = cerByParN(goodGenes);
y = parByCerN(goodGenes);

subplot(1,3,3);
scatter(x, y, sz, [.7 .7 .7]);
hold on;
scatter(x(gList), y(gList), sz, z(gList), 'filled', 'MarkerEdgeColor', 'k');
xlabel('Rc, hyb.cer - NNhyb.par');
ylabel('Rc, hyb.par - NNhyb.cer');
axis square; grid on; box on; set(gca, 'fontsize', fs);
text(x(exampleGenesIdx), y(exampleGenesIdx), listOfGenes(exampleGenesIdx), 'fontsize', fs);
axis([-1 1 -1 1]);
set(gca, 'xtick', -1:.5:1, 'ytick', -1:.5:1);
cb = colorbar; cb.Label.String = 'Rc, hybrid alleles'; cb.FontSize = fs;
set(gcf, 'color', 'w');
% text(x(gList), y(gList), listOfGenes(gList));

%% cer-like, par-like
x = obs_all_gg(:, 3); % cer-hyc
y = obs_all_gg(:, 4); % par-hyp
cerLike1 = find(x(gList) - y(gList) > diffTh);
parLike1 = find(y(gList) - x(gList) > diffTh);
additive = find(x(gList) - y(gList) < diffTh & x(gList) - y(gList) > -diffTh);
genes = [cerLike1; parLike1; additive];
genesNameIdx = intint(goodGenes(gList(genes)));
geneNames = allGenes(genesNameIdx);
category = vertcat(repmat({'cer-like'}, length(cerLike1), 1), repmat({'par-like'}, length(parLike1), 1), ...
    repmat({'additive'}, length(additive), 1));
Rc_scores = obs_all_gg(gList(genes), :);
NNwithinSpecies = nnWithinGenome.secondOrder.corrs(goodGenes(gList(genes)), 1:2);
NNcrossSpecies = NNtop5.corr(goodGenes(gList(genes)), 1:2);
t = table(geneNames, category, Rc_scores(:, 1), Rc_scores(:, 2), Rc_scores(:, 3), Rc_scores(:, 4), ...
    NNwithinSpecies(:, 1), NNwithinSpecies(:, 2), NNcrossSpecies(:, 1), NNcrossSpecies(:, 2), ...
    GP.gene_infoR64.desc(genesNameIdx), ...
    'variableNames', {'GeneName', 'category', 'Rc_cer-par', 'Rc_hybrid', 'Rc_cer-hybrid_par', 'Rc_par-hybrid_par', ...
    'Rc_cer_NNcer', 'Rc_par_NNpar', 'Rc_cer_NNpar', 'Rc_par_NNcer', ...
    'SGD_description'});
t = sortrows(t, [2, 3]);
writetable(t, [homeDir, 'revision/cis_cerLike_parLike.xlsx']);

disp(['cer like : ', num2str(length(cerLike1)), ', ', num2str(length(cerLike1)*100./length(gList)), ' %']);
disp(['par like : ', num2str(length(parLike1)), ', ', num2str(length(parLike1)*100./length(gList)), ' %']);
disp(['additive : ', num2str(length(additive)), ', ', num2str(length(additive)*100./length(gList)), ' %']);