function go_process_out(frootname) % Reformat results from FUNC GO analysis. % % Input: % frootname: GO result root name % % go_batch.m % Written by Andreas Wollstein, Ludwig Maximilians University of Munich. % Copyright (c) 2015 % % This file is part of a set of Matlab scripts published with: % Bozicevic et al. 2015 Molecular Ecology % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program. If not, see . f = dir(['go_FUNC/' frootname '/refine*']); CC = {}; if length(f)==3 [cc1 descr1 go1 sig1 pu1 po1 pru1 pro1 xx1] = textread(['go_FUNC/' frootname '/' f(1).name],'%s %s %s %s %s %s %s %s %s','delimiter','\t','headerlines',1); [cc2 descr2 go2 sig2 pu2 po2 pru2 pro2 xx1] = textread(['go_FUNC/' frootname '/' f(2).name],'%s %s %s %s %s %s %s %s %s','delimiter','\t','headerlines',1); [cc3 descr3 go3 sig3 pu3 po3 pru3 pro3 xx1] = textread(['go_FUNC/' frootname '/' f(3).name],'%s %s %s %s %s %s %s %s %s','delimiter','\t','headerlines',1); ix = strcmp(sig1,'-'); cc1(ix) = []; descr1(ix) = []; go1(ix) = []; pru1(ix) = []; pro1(ix) = []; ix = strcmp(sig2,'-'); cc2(ix) = []; descr2(ix) = []; go2(ix) = []; pru2(ix) = []; pro2(ix) = []; ix = strcmp(sig3,'-'); cc3(ix) = []; descr3(ix) = []; go3(ix) = []; pru3(ix) = []; pro3(ix) = []; res.cc = [cc1; cc2 ;cc3]; res.descr = [descr1 ; descr2; descr3]; res.go = [go1 ; go2 ; go3]; res.goid = str2double(strrep(res.go,'GO:','')); res.pru = str2double([pru1 ; pru2; pru3]); res.pro = str2double([pro1 ; pro2; pro3]); res.pu = str2double([pu1 ; pu2; pu3]); res.po = str2double([po1 ; po2; po3]); mp = min(res.pro,res.pru); [s idx] = sort(mp); res.cc = res.cc(idx); res.descr = res.descr(idx); res.go = res.go(idx); res.goid = res.goid(idx); res.pru = res.pru(idx); res.pro = res.pro(idx); res.pu = res.pu(idx); res.po = res.po(idx); save(['go_FUNC/' frootname '.mat'],'res'); %% save top ten f = fopen(['go_FUNC/' frootname '.t10.txt'],'w'); fprintf(f,'Category\tName\tId\tP_u\tP_o\n'); m = min(10,length(res.pro)); for i=1:m fprintf(f,'%s\t%s\t%s\t%e\t%e\n',res.cc{i},res.descr{i},res.go{i},res.pru(i),res.pro(i)); end fclose(f); %% save top twenty f = fopen(['go_FUNC/' frootname '.t20.txt'],'w'); fprintf(f,'Category\tName\tId\tP_u\tP_o\n'); m = min(20,length(res.pro)); for i=1:m fprintf(f,'%s\t%s\t%s\t%e\t%e\n',res.cc{i},res.descr{i},res.go{i},res.pru(i),res.pro(i)); end fclose(f); %% save bonf < 0.01 mp = min([res.pru,res.pro],[],2) .* 493076; idx = find(mp<0.01); f = fopen(['go_FUNC/' frootname '.bonf.txt'],'w'); fprintf(f,'Category\tName\tId\tP_u\tP_o\tP_bonf\n'); for i=1:length(idx) fprintf(f,'%s\t%s\t%s\t%e\t%e\t%e\n',res.cc{idx(i)},res.descr{idx(i)},res.go{idx(i)},res.pru(idx(i)),res.pro(idx(i)),mp(idx(i))); end fclose(f); %% save all f = fopen(['go_FUNC/' frootname '.all.txt'],'w'); fprintf(f,'Category\tName\tId\tP_u\tP_o\n'); for i=1:length(res.pro) fprintf(f,'%s\t%s\t%s\t%e\t%e\n',res.cc{i},res.descr{i},res.go{i},res.pru(i),res.pro(i)); end fclose(f); %% output generic_format for enrichment map %% ID DESCRIPTION P-VALUE FDR PHENOTYPE f = fopen(['go_FUNC/' frootname '.enrmap.txt'],'w'); mp = min([res.pru,res.pro],[],2); fprintf(f,'ID\tDESCRIPTION\tP-VALUE\tFDR\tPHENOTYPE\n'); idx = find(mp<0.01); for i=1:length(idx); pheno=1; if res.pru(idx(i))<0.01 pheno=-1; end fdr = min(res.pru(idx(i)),res.pro(idx(i))); fprintf(f,'GO:%07d\t%s\t%f\t%f\t%d\n',res.goid(idx(i)),res.descr{idx(i)},fdr,fdr,pheno); end fclose(f); fprintf('all done!\n'); %% output generic_format for enrichment map, BONFERRONI Corrected %% ID DESCRIPTION P-VALUE FDR PHENOTYPE mp = min([res.pru,res.pro],[],2) .* 493076; f = fopen(['go_FUNC/' frootname '.enrmap.bonf.txt'],'w'); fprintf(f,'ID\tDESCRIPTION\tP-VALUE\tFDR\tPHENOTYPE\n'); idx = find(mp<0.01); for i=1:length(idx) pheno=1; if res.pru(idx(i))<0.01 pheno=-1; end fdr = min(res.pru(i),res.pro(idx(i))); fprintf(f,'GO:%07d\t%s\t%f\t%f\t%d\n',res.goid(idx(i)),res.descr{idx(i)},fdr,fdr,pheno); end fclose(f); fprintf('all done!\n'); else fprintf('%s didnt work!\n',frootname); end end