clear;

%% Pi histogram
load p.txt;

figure; set(gcf,'Position',[793,453,356,288]);
[n,x] = hist(p,30); f=n/sum(n*(x(2)-x(1))); bar(x,f,'histc'); 
box off; colormap([[1 1 1];[0 0 0]]);
xlim([0,1.02]);
set(gca,'Position',[0.13,0.13,0.775,0.795],'FontSize',14);
set(get(gca,'Children'),'LineWidth',1.15);

xlabel('$$\pi_g$$','FontSize',16,'Interpreter','latex');
ylabel('Frequency','FontSize',14);
box on;

%% Report top pathways
% load information
fid = fopen('geneset_id.txt','r');
x = textscan(fid,'%s','delimiter','\n','BufSize',50000);
fclose(fid);
x = x{1};
J = length(x);
for j = 1:J
  g = strread(x{j});
  geneset_id{j} = g(1:end-1)';
  geneset_length(j) = length(geneset_id{j});
end
geneset_id = geneset_id';
clear x g;

geneset_name = textread('geneset_name.txt','%s','delimiter','\n');

load geneid.txt;

load p_sign.txt;

load result.txt;
LB = result(:,1); UB = result(:,2);

% plot Log ML bounds (sort by UB)
figure; set(gcf,'Position',[793,453,356,288]);
[tmp,ind] = sort(UB,'descend');
hL = plot(LB(ind),'x','MarkerEdgeColor','k','MarkerSize',6);
hold on;
hU = plot(tmp,'o','MarkerEdgeColor','k','MarkerSize',6);
xlim([-50,960]); 
xlabel('$$j$$','FontSize',16,'Interpreter','latex');
ylabel('Log marginal likelihood','FontSize',14);
set(gca,'Position',[0.18,0.13,0.775,0.795],'FontSize',14);
set(get(gca,'Children'),'LineWidth',1.15);
h = legend([hU,hL],'Upper bound','Lower bound');
legend boxoff;
set(h,'Location','NorthEast','FontSize',14);

plot(LB(ind),'k:.','MarkerEdgeColor','k','MarkerSize',1);
plot(tmp,'k:.','MarkerEdgeColor','k','MarkerSize',1);

% select top pathways and generate report
top = 30;
nump = length(UB);
[sortlogml,index] = sort(UB);
sortlogmlLB = LB(index);

toppathways = [];
for i = 1:top
  j = nump-i+1;
  ind = index(j);
  [c,ia,ib] = intersect(geneset_id{ind},geneid);
  
  pwj = struct('Rank',i,'Name',geneset_name{ind},'Description','',...
    'LogMLupperB',sortlogml(j),'LogMLlowerB',sortlogmlLB(j),...
    'GeneSetSize',geneset_length(ind),'Reliability',sum(p(ib)>min(p))>0,...
    'GeneSetEntrezId',c,...
    'AveAssoProb',p(ib)');
  toppathways = [toppathways;pwj];
end

fname = ['top',int2str(top),'genesets_ER.txt'];
fid = fopen(fname,'w');
fprintf(fid,'%s\t%s\t%s\t%s\t%s\t%s\t%s\n',...
  'Index', 'Pathway','Size',...
  'Log ML (UB)','Log ML (LB)','Bound Dist.','');
  
for i = 1:top
  j = nump-i+1;
  trust = '';
  if toppathways(i).Reliability == 0
    trust = '*';
  end
  fprintf(fid,'%d\t%s\t%d\t%5.2f\t%5.2f\t%5.2f\t%s\n', ...
    toppathways(i).Rank,toppathways(i).Name, ...
    toppathways(i).GeneSetSize,...
    toppathways(i).LogMLupperB,toppathways(i).LogMLlowerB,...
    toppathways(i).LogMLupperB-toppathways(i).LogMLlowerB,...
    trust);
end
fprintf(fid,'\n(The pathways marked with * are not reliable)');
fclose(fid);

%% Plot interesting pathways
jid = [1,2,4,6,11,12,14,18,29]; J = length(jid);
xltxt = 'abcdefghijkl';
figure; set(gcf,'Color','w','Position',[292,178,754,510]);
for j = 1:J
  subplot(3,3,j);
  [tmp,ia,ib] = intersect(toppathways(jid(j)).GeneSetEntrezId,geneid);
  hist(p(ib(find(p_sign(ib) == -1))),30);
  hneg = findobj(gca,'Type','patch');
  set(hneg,'FaceColor',[1,1,1]);
  xlim([0,1]);
  hold on;
  hist(p(ib(find(p_sign(ib) == 1))),30);
  h = findobj(gca,'Type','patch');
  hpos = setxor(hneg,h);
  set(hpos,'FaceColor',[0.6,0.6,0.6],'FaceAlpha',0.7);
  xlabel(['(',xltxt(j),')'],'FontSize',12);
  pos = get(gca,'Position');
  pos(1) = 0.04+mod(j-1,3)/3; pos(3) = 0.25;
  pos(2) = 0.75-0.32*floor((j-1)/3); pos(4) = 0.21;
  set(gca,'Position',pos);
  pname  = toppathways(jid(j)).Name;
  ind = findstr('_',pname);
  pname(ind) = ' ';
  title(['A',int2str(jid(j)),': ',pname]);%,toppathways(j).Name]);
  hl = legend([hpos,hneg],'Positive','Negative');
  legend boxoff;
  set(hl,'Location','NorthWest','FontSize',12);
end