% figure S6
%% A
% CHA1 in external datasets
% load Gasch and Michal
load([homeDir, 'Ext_data/Gasch/gasch.mat']);
run([homeDir, 'gene_gene_corr/old scripts/Michal_MSN_order.m']);
%%
lowN_michal = find(contains(Michal.Labels, 'LowNitrogen'));
lowN_gasch = find(contains(gasch.conds, 'Nitrogen Depletion'));
aa_gasch = find(contains(gasch.conds, 'Amino acid'));

%%
currGene = 'CHA1';
gn = GP.gene_table.(currGene);

dGasch = gasch.data_id(gn, :);
dMichal = Michal.Data(gn, :);
lowN = {lowN_aa_gasch, lowN_michal};
ylabs = {'log2 fold change', 'log2 expression'};
figure;
currD = dGasch;
[~, sortdidx] = sort(currD);
subplot(1,2,1); hold on;
scatter(1:length(sortdidx), currD(sortdidx));
lowN_idx = ismember_smart(sortdidx, lowN_gasch);
scatter(lowN_idx, currD(sortdidx(lowN_idx)));
aa_idx = ismember_smart(sortdidx, aa_gasch);
scatter(aa_idx, currD(sortdidx(aa_idx)));
legend({'Other conditions', 'Nitrogen Depletion', 'Amino Acid + adenine starvation'});
xlabel('sorted samples');
ylabel(ylabs{1});
axis square; box on; grid on; set(gca, 'fontsize', fs);

subplot(1,2,2); hold on;
currD = dMichal;
[~, sortdidx] = sort(currD);
scatter(1:length(sortdidx), currD(sortdidx));
lowN_idx = ismember_smart(sortdidx, lowN_michal);
scatter(lowN_idx, currD(sortdidx(lowN_idx)));
legend({'Other conditions', 'Low Nitrogen'});
xlabel('sorted samples');
ylabel(ylabs{2});
axis square; box on; grid on; set(gca, 'fontsize', fs);

set(gcf, 'color', 'w');    
%% B 
% currGene = 'CHA1';
% gn = GP.gene_table.(currGene);
GCwt = 404:418;
GChap4 = 419:433;
time_after_start_h = [0    1.3333    2.0833    2.8333    3.5833    4.3333    5.0833    6.0833    7.0833    8.0833    9.0833   10.0833   11.0833   23.1667   32.5000];
figure; 
hold on;
colors = lines(4);
h = cell(4,1);
for i = 1:4
%     d = [expStruct.(gb{i}).exp(gn, [GCwt]); expStruct.(gb{i}).exp(gn, [GChap4])];
    d = expStruct.(gb{i}).exp(gn, [GCwt]);
    h{i} = plot(time_after_start_h(1:13), d(1:13), 'o', 'color', colors(i, :));
    plot(time_after_start_h(1:13), nanmean(d(1:13), 1), 'color', colors(i, :));
end
xlabel('time (hours)');
ylabel('log2 expression');
legend([h{1}, h{2}, h{3}, h{4}], {'cer', 'par', 'hyb.cer', 'hyb.par'}, 'location', 'southeast');
axis square; box on; grid on; set(gca, 'fontsize', fs);
set(gcf, 'color', 'w');
%% Michals entry exit stationary
figure;
d = MSN.SCS_30c.expression(gn, :, 1);
dtime = MSN.SCS_30c.time;
% dtime = [dtime(1:21), dtime(22:end)-70];
subplot(1,2,1);
plot(dtime(1:20), d(1:20), '.-');
xlabel('time (hours)');
ylabel('log2 expression');
 box on; grid on; set(gca, 'fontsize', fs);
 ylim([8.5 11]);
 
subplot(1,2,2);
plot(dtime(22:end), d(22:end), '.-');
xlabel('time (hours)');
ylabel('log2 expression');
 ylim([8.5 11]);
 box on; grid on; set(gca, 'fontsize', fs);
set(gcf, 'color', 'w');


%% cha4 expression in WT
gn = GP.gene_table.CHA4;
expWT = cell(3,4);
for i = 1:length(TPs)
    currTP = TPs{i};
    for j = 1:length(gb)
        currEntry = ['WCD', gb{j}];
        idxSamples = find(contains(dataLabels.(currEntry), ['WT ', currTP]));
        expWT{i, j} = dataS.(currEntry)(gn, idxSamples);
    end
end

figure; 
for i = 1:3
    subplot(1,3,i); hold on;
    means = cellfun(@mean, expWT(i, :));
    stds = cellfun(@std, expWT(i, :));
    bar(means, 'FaceColor', rgb('lightblue'));
    errorbar(means, stds, '.b');
    plotSpread(expWT(i, :));
    [h, p, ci, stats] = ttest2(expWT{i, 1}, expWT{i, 2});
    disp(['TP', num2str(i)]);
    disp(['cer-par : p-value ', num2str(p)]);
    disp(['t-stat: ', num2str(stats.tstat)]);
    disp(['effect size: ', num2str(nanmean(expWT{i, 1}) - nanmean(expWT{i, 2}))]);
    set(gca, 'xtick', 1:4, 'xticklabels', {'cer', 'par', 'hyb.cer', 'hyb.par'}); xtickangle(90);
    box on; 
    ylim([0 7]);
    ylabel('log2 expression');
    set(gca, 'fontsize', fs);
end
set(gcf, 'color', 'w');

