function results = GOTermFinder(genes,background_genes,pvalCutoff,ontology),
% GOTermFinder    Uses hypergeometric distribution of gene group given a
%                 background set of genes.  P-values are
%                 Bonferroni-corrected to account for multiple hypotheses.
%
% Syntax:
%   results=GOTermFinder(genes,background_genes,pvalCutoff,ontology)
%
% Inputs:
%   genes - cell array of gene group to be tested (must be IDs corresponding present all_commonNames 
%            vector)
%   background_genes - background gene group, (e.g. genes present on array)
%   pvalCutoff - significance level (suggested .05)
%   ontology - GO_process, GO_function, GO_component as contained in
%               GO_matlab or generated by create_GO_reference
%
% Outputs:
%   results -
%
% Author: Chad L. Myers (clmyers@princeton.edu)
% Date: 11/7/05


    genes = unique(genes);
    background_genes = unique(background_genes);
    
    gene_orfs = ontology.all_commonNames;
    set_ids=zeros(length(genes),1);
    
    geneIDHash = java.util.HashMap;
    for i=1:length(ontology.all_commonNames),
        geneIDHash.put(ontology.all_commonNames{i},i);
    end
    
    for i=1:length(genes),
       if(geneIDHash.containsKey(genes{i}))
           set_ids(i) = geneIDHash.get(genes{i});
       else 
           set_ids(i) =-1;
           fprintf('Warning! Could not recognize test gene: %s!\n',genes{i});
       end
    end
    
    ind = find(set_ids > 0);
    set_ids=set_ids(ind);    

 if(length(background_genes) == 0) background_ids = 1:length(ontology.all_commonNames);
 else    
    background_ids=zeros(length(background_genes),1);
    for i=1:length(background_genes),
      if(geneIDHash.containsKey(background_genes{i}))
           background_ids(i) = geneIDHash.get(background_genes{i});
       else 
           background_ids(i) =-1;
           %fprintf('Warning! Could not recognize background gene: %s!\n',background_genes{i});
       end
    end
 end
    ind = find(background_ids > 0);
    background_ids=background_ids(ind);    
    
    results=cell(1);
    pvals=[];
    
    ind=find(sum(ontology.annotation_matrix(set_ids,:),1) > 1);
    allNodes=ind;
    
    total_genes = length(background_genes);
    %total_genes = length(background_ids);
    
    allNodes = unique(allNodes);
    count=1;
    for i=1:length(allNodes),
        %pval = 1-hygecdf(sum(ontology.annotation_matrix(set_ids,allNodes(i)))-1,total_genes,sum(ontology.annotation_matrix(background_ids,allNodes(i)))+0,length(set_ids));
        pval = 1-hygecdf(sum(ontology.annotation_matrix(set_ids,allNodes(i)))-1,total_genes,sum(ontology.annotation_matrix(background_ids,allNodes(i)))+0,length(genes));

        if(pval < pvalCutoff/length(allNodes))
            %get list of genes belonging to this node
            ind = find(ontology.annotation_matrix(set_ids,allNodes(i)) > 0);
            temp =[];
            for j=1:length(ind),
                temp = [temp,ontology.all_commonNames{set_ids(ind(j))}];
                if(j <length(ind)) temp =[temp,', '];end
            end
            results{count}= sprintf('%s\t%s\t%d of %d (%.0f%%)\t%d of %d (%.1f%%)\t%.3e\t%s',ontology.all_nodes{allNodes(i)},ontology.all_nodes_desc{allNodes(i)},sum(ontology.annotation_matrix(set_ids,allNodes(i)))+0,length(genes),100*(sum(ontology.annotation_matrix(set_ids,allNodes(i)))+0)/length(genes),(sum(ontology.annotation_matrix(background_ids,allNodes(i)))+0),total_genes,100*(sum(ontology.annotation_matrix(background_ids,allNodes(i)))+0)/total_genes,pval*length(allNodes),temp);
            %results{count}= sprintf('%s\t%s\t%d of %d (%.0f%%)\t%d of %d (%.1f%%)\t%.3e\t%s',ontology.all_nodes{allNodes(i)},ontology.all_nodes_desc{allNodes(i)},sum(ontology.annotation_matrix(set_ids,allNodes(i)))+0,length(set_ids),100*(sum(ontology.annotation_matrix(set_ids,allNodes(i)))+0)/length(set_ids),(sum(ontology.annotation_matrix(background_ids,allNodes(i)))+0),total_genes,100*(sum(ontology.annotation_matrix(background_ids,allNodes(i)))+0)/total_genes,pval*length(allNodes),temp);
            pvals(count)= pval;
            count=count+1;    
        end
    end

[vals,ind]=sort(pvals);
results=results(ind);

%print GO term results
for i=1:length(results),
    fprintf('%s\n',results{i});
end
