%% how many genes change per condition?

DE_sp = load([homeDir, 'revision/DESeq_normalization/DE.mat'], 'spGenes');
DE_trans = load([homeDir, 'revision/DESeq_normalization/DE.mat'], 'transGenes');
DE_cis = load([homeDir, 'revision/DESeq_normalization/DE.mat'], 'cisGenes');
DE_sp = DE_sp.spGenes;
DE_trans = DE_trans.transGenes;
DE_cis = DE_cis.cisGenes;

currLfcTh = 1;

genes_cond_sp = abs(comparable_samples_de.sp) > currLfcTh;
genes_cond_cis = abs(comparable_samples_de.cis) > currLfcTh;
genes_cond_trans = abs(comparable_samples_de.trans) > currLfcTh;

% genes_per_cond = table(sum(genes_cond_sp, 1)', sum(genes_cond_cis, 1)', sum(genes_cond_trans, 1)', 'variableNames', {'sp', 'cis', 'trans'});
genes_per_cond = table(sum(genes_cond_sp(DE_sp, :), 1)', sum(genes_cond_cis(DE_cis, :), 1)', ...
    sum(genes_cond_trans(DE_trans, :), 1)', 'variableNames', {'sp', 'cis', 'trans'});
bw=50;
figure; histogram(genes_per_cond.sp, 'binwidth', bw); hold on;
histogram(genes_per_cond.cis, 'binwidth', bw);
histogram(genes_per_cond.trans, 'binwidth', bw);
xlabel('# differentially-expressed genes')
ylabel('# comparable samples');
legend({'overall', 'cis', 'trans'});
fs = 14;
set(gca, 'fontsize', fs);
set(gcf, 'color', 'w');
% %% In how many conditions does a gene change?
% 
% cond_per_gene_mat = nan(6701, 3);
% cond_per_gene_mat(DE_sp, 1) = sum(genes_cond_sp(DE_sp, :), 2);
% cond_per_gene_mat(DE_cis, 2) = sum(genes_cond_cis(DE_cis, :), 2);
% cond_per_gene_mat(DE_trans, 3) = sum(genes_cond_trans(DE_trans, :), 2);
% cond_per_gene = table(cond_per_gene_mat(:, 1), cond_per_gene_mat(:, 2), cond_per_gene_mat(:, 3), ...
%     'variableNames', {'sp', 'cis', 'trans'});
% bw=50;
% figure; histogram(cond_per_gene.sp, 'binwidth', bw); hold on;
% histogram(cond_per_gene.cis, 'binwidth', bw);
% histogram(cond_per_gene.trans, 'binwidth', bw);
% xlabel('# conditions'); ylabel('# genes');

%%
% effects = {'sp', 'cis', 'trans'};
% figure; hold on;
% for i = 1:length(effects)
%     d = cond_per_gene.(effects{i});
%     [ycdf, xcdf] = cdfcalc(d);
%     yccdf = 1 - ycdf(1:end-1);
% %     yccdf = yccdf*length(find(~isnan(d)));
%     plot(xcdf, yccdf)
% end
% legend(effects)
%% example genes
cis_trans_specific_examples3;
%% strongest effects
N= 200;
sumCis = sum(abs(comparable_samples_de.cis(DE_cis, :)), 2);
sumTrans = sum(abs(comparable_samples_de.trans(DE_trans, :)), 2);
[~, sortdCis] = sort(sumCis, 'descend');
topCis = sortdCis(1:N);
[~, sortdTrans] = sort(sumTrans, 'descend');
topTrans = sortdTrans(1:N);

%
wantedLists = struct;
% wantedLists.cis = DE_cis;
% wantedLists.trans = DE_trans;
wantedLists.cis = DE_cis(topCis);
wantedLists.trans = DE_trans(topTrans);
f = {'cis', 'trans'};
figure;
% locs = {1:9, 10, 13:21, 22};
% sprange = [1 22];
locs = {1 2};
sprange = [1 2];
n=1;
cm1 = cbrewer2('RdBu'); cm1 = cm1(end:-1:1, :);
cm2 = cbrewer2('Purples');
orders = cell(1,2);
cax = [-5, 5];
for i = 1:2
    currList = wantedLists.(f{i});
    currDraw = comparable_samples_de.(f{i})(currList, :);
%     pd = pdist(currDraw);
%     l = linkage(pd);
%     sortdidx = optimalleaforder(l, pd);
    cgo = clustergram(currDraw, 'colormap', redbluecmap, 'DisplayRange', 5);
    gorder = cellfun(@(x) str2num(x), cgo.RowLabels);
    corder = cellfun(@(x) str2num(x), cgo.ColumnLabels);
    orders{i} = gorder;
    clear cgo;
    
    currDraw = currDraw(gorder, corder);
    s1 = subplot(sprange(1), sprange(2), locs{n});
    n = n+1;
    imagesc(currDraw);
    colormap(s1, cm1); caxis(cax);
    xlabel('Conditions'); ylabel('Genes');

%     s2 = subplot(sprange(1), sprange(2), locs{n});
%     flex = mad(currDraw, [], 2);
%     n = n+1;
%     imagesc(flex);
%     colormap(s2, cm2);
%     caxis([0 10]);
%     set(gca, 'yticklabels', {}); set(gca, 'xticklabels', {}); 
end
set(gcf, 'color', 'w', 'position', [474   204   784   307]);
%% both
both = union(wantedLists.cis, wantedLists.trans);
toCluster = comparable_samples_de.trans(both, :);
pd = pdist(toCluster);
l = linkage(pd);
sortdidx = optimalleaforder(l, pd);

figure;
for i = 1:2
    toDraw = comparable_samples_de.(f{i})(both(sortdidx), :);
    subplot(1,2,i);
    imagesc(toDraw);
    colormap(cm1); caxis(cax);
    xlabel('Conditions'); ylabel('Genes');
end

