GO.h

/**************************************************************************
 * Copyright (c) 2004-2011 T. M. Murali                                   *
 * Copyright (c) 2007-2011 Naveed Massjouni                               *
 * Copyright (c) 2009-2011 Christopher L. Poirel                          *
 * Copyright (c) 2011 Phillip Whisenhunt                                  *
 * Copyright (c) 2011 David Badger                                        *
 * Copyright (c) 2010 Clifford Conley Owens III                           *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm 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.                                    *
 *                                                                        *
 * Biorithm 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 Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/


#ifndef _GO_H
#define _GO_H


#include<iostream>
#include<map>
#include<set>
#include<vector>
#include<string>
#include<sstream>

#include "boost/tuple/tuple_io.hpp"
#include "biofunction.h"
#include "enrichment.h"
// HACKHACKHACK. this #include is just for defining the MyNT type.
#include "histogram.h"
#include "setTemplates.C"

using namespace std;

// forward declaration
class MyAnnotations;

//enumated type for GO relationships; 0-is_a_or_part_of, 1-is_a, 2-part_of, 3-regulates, 4-negatively_regulates, 5-positively_regulates
enum enumGoRelationships
{
    is_a_or_part_of,
    is_a,
    part_of,
    regulates,
    negatively_regulates,
    positively_regulates
};

class GOFunction  : public BioFunction {
  friend class GeneOntology;
        
public:
  GOFunction()
      : BioFunction()
    {
      _init();      
    }

  GOFunction(string cat, string id, string name = "")
      : BioFunction(cat, id, name)
    {
      _init();
    }

  virtual ~GOFunction()
    {}
  

  // set up all the internal variables.
  void _init()
    {
      integerId=-1;
//      name=name_space=
      def=comment="";
      obsolete=false;
      _siblingsComputed[0] =_siblingsComputed[1] =
        _siblingsComputed[2] = false;
      
    }

  MyNT computeSemanticValue() const;
  
  
  unsigned int computeSiblings(unsigned int type = 0);
  



  int getIntegerId() {return integerId;}

  vector<int> getAltIds() {return alt_id;}

  vector<string> getSubsets() {return subset;}

  vector<string> getSynonyms() {return synonym;}

  vector<string> getXrefAnalogs() {return xref_analog;}

//   /*! Accessor for name
//    * \returns function name*/
//   string getName() {return name;}

  string getNamespace()
    {return getCategory();}
//    {return name_space;}

  string getDef() {return def;}

  string getComment() {return comment;}

  bool isObsolete() {return obsolete;}
  
  virtual unsigned int getNumParents(int type = 0)
    {
      return(parents[type].size());
    }
  virtual unsigned int getNumChildren(int type = 0)
    {
      return(children[type].size());
    }
  
  virtual set< GOFunction *> getAncestors(int type = 0)
    {
      
      //check if already computed/no parents
      if(ancestors[type].size()>0)
        return ancestors[type];
      ancestors[type]= parents[type];
      for(set< GOFunction *>::iterator i = parents[type].begin();
          i != parents[type].end(); i++)
        ancestors[type]= setunion((*i)->getAncestors(type), ancestors[type]);
  
      return ancestors[type];
    }

  virtual map< GOFunction *, unsigned int > getAncestorsWithDistance(int type = 0);
  
  virtual set< GOFunction *> getChildren(int type = 0)
    {
      return(children[type]);
    }
  
  virtual set< GOFunction *> getDescendants(int type = 0)
    {
      //check if already computed or no descendants.
      if(descendants[type].size()>0)
        return descendants[type];
      descendants[type] = children[type];
      for(set< GOFunction *>::iterator i = children[type].begin();
          i != children[type].end(); i++)
        descendants[type]= setunion((*i)->getDescendants(type), descendants[type]);
  
      return descendants[type];
    }
  
  virtual map< GOFunction *, unsigned int > getDescendantsWithDistance(int type = 0);

  virtual unsigned int getDistanceToAncestor(GOFunction *ancestor,
                                             unsigned int type = 0);
  
  
  virtual unsigned int getDistanceToDescendant(GOFunction *descendant,
                                               unsigned int type = 0);
  
  
  virtual set< GOFunction *> getParents(int type = 0)
    {
      return(parents[type]);
    }
  
  set< GOFunction * > getSiblings(unsigned int type = 0)
    {
      if (!_siblingsComputed[type])
        computeSiblings(type);
      return(_siblings[type]);
    }
  

  bool isAncestor(GOFunction *candidate) 
    {
      set< GOFunction * > ancestors = getAncestors();
      return(ancestors.end() != ancestors.find(candidate));
    }

  bool isChild(GOFunction* candidate)
    {
      set< GOFunction* > children = getChildren();
      return(children.end() != children.find(candidate));
    }

  bool isDescendant(GOFunction *candidate) 
    {
      set< GOFunction * > descendants = getDescendants();
      return(descendants.end() != descendants.find(candidate));
    }

  bool isParent(GOFunction* candidate)
    {
      set< GOFunction* > parents = getParents();
      return(parents.end() != parents.find(candidate));
    }

  
  void addDepth(unsigned int depth, int type = 0)
    {
      _depths[type].insert(depth);
    }

  set< unsigned int > getDepths(int type = 0)
    {
      return(_depths[type]);
    }

  unsigned int getMinimumDepth(int type = 0)
    {
      return(*(_depths[type].begin()));
    }
  
  unsigned int getMaximumDepth(int type = 0)
    {
      set< unsigned int >::iterator ditr = _depths[type].end();
      ditr--;
      return(*ditr);
    }
  

protected:
  unsigned int _computeSiblingsNaively(unsigned int type = 0);

  
  
private:
  int integerId;
  vector<int> alt_id;
  vector<string> subset,synonym,xref_analog;
  string def,comment;
  bool obsolete;

  //TODO: Make these enumerated instead of 0,1,2
  vector<int> tempParents[6];                   //0-all, 1-is_a, 2-part_of
  vector<int> tempChildren[6];                  //0-all, 1-is_a, 2-part_of
  set<GOFunction*> parents[6],ancestors[6];     //0-all, 1-is_a, 2-part_of
  set<GOFunction*> children[6], descendants[6]; //0-all, 1-is_a, 2-part_of
  set<GOFunction*> _siblings[6];
  bool _siblingsComputed[6];
  map< GOFunction*, unsigned int > ancestorsWithDistance[6], descendantsWithDistance[6];

  
  // other computed info.
  set< unsigned int > _depths[3];
};

enum GOGroupByDepthType {
  GO_GROUP_BY_MINIMUM_DEPTH = 0, GO_GROUP_BY_AVERAGE_DEPTH,
  GO_GROUP_BY_MAXIMUM_DEPTH 
};

class GeneOntology {

public:
  GeneOntology()
    : default_namespace(""), functions(), idxId(), idxName(), _functionCategoryAbbrevs(),
      _legalFunctionCategories(), _knownFunctionIds(), _topologicalSortComputed(false),
      _topologicallySortedFunctions()
    {
      _init();
    }

  GeneOntology(istream *in);
  
  void read(string file);
  void read(istream *in);

  void readTerms(string file, vector<GOFunction *> &functions);

  void computeEnrichmentsOntologizer(const set< string > &geneSet,
                                     // not const becaus MyAnnotations::numAnnotatedGenes() is not a const function.
                                     MyAnnotations &annotations,
                                     vector<EnrichmentRecord< string, string> > &rejalts) const;
  

  void computeFunctionsWithNoChildren(set< string > &functions,
                                      set< string > &childless);

  void computeFunctionsWithNoChildren(set< GOFunction* > &functions,
                                      set< GOFunction* > &childless);

    

  typedef boost::tuple< GOFunction *, unsigned int, unsigned int > GOLCATuple;
  int computeLeastCommonAncestors(GOFunction *function1, GOFunction *function2,
                                  vector< GOLCATuple >&lcas);
//                                  vector< GOFunction * > &lcas);


  void computeMostSpecificFunctions(set< GOFunction * > &functions,
                                    set< GOFunction * > &specificFunctions);
  

  template < typename InputIterator, typename OutputIterator, typename Extractor >
  void computeMostSpecificFunctionsTemplated(
    InputIterator functionsBegin, InputIterator functionsEnd,
    OutputIterator specificFunctionsBegin,
    Extractor _extractor
    );

  void computeMostSpecificFunctions(set< string > &functions,
                                    set< string > &specificFunctions);
  


  set<GOFunction*> getAncestors(GOFunction* obj, int type=0);
  
  /* Setter for boolRelationship
   * \param r used to set boolRelationship */
  void setBoolRelationship(bool r){_boolRelationship = r;}
  
  /* Setter for boolRelationship
   * \param r used to set boolRelationship */
  void setStrRelationshipFile(string file){_strRelationshipFile = file;}
  
  map< GOFunction*, unsigned int >getAncestorsWithDistance(GOFunction* obj,
                                                           int type = 0);

  set<GOFunction*> getParents(GOFunction* obj, int type=0);

  set<GOFunction*> getDescendants(GOFunction* obj, int type=0);


  map< GOFunction*, unsigned int > getDescendantsWithDistance(GOFunction* obj, int type = 0);
  
  
  set<GOFunction*> getChildren(GOFunction* obj, int type=0);
  
  GOFunction* getFunctionById(int id) const;
  
  GOFunction* getFunctionById(string id) const
    {
      unsigned int intID;
      stringstream ss;
      ss << id;
      ss >> intID;
      return(getFunctionById(intID));
    }

  GOFunction* getFunctionByName(string name) const;

  /* Return true if the given id is known by the instance of GO.
     Otherwise, return false.
   */
  bool isKnownFunctionId(string id) const
    {
      if (_knownFunctionIds.find(id)==_knownFunctionIds.end())
        return false;
      else
        return true;
    }

  
  void computeFunctionDepths(//vector< set< GOFunction* > > &functionsByDepth,
    unsigned int type = 0);

  void computeSiblings(unsigned int type = 0);
  
  void computeTopologicalSort(vector< GOFunction* > &sortedFunctions, unsigned int type = 0);
  

  void printFunctionDepths(ostream &fstr, unsigned int type = 0);

  void printFunctionDescendants(ostream &fstr, unsigned int type = 0);

  void printRelationships();
  
  void groupFunctionsByDepth(map< unsigned int, set< GOFunction * > > &functionsByDepth,
                             string functionCategory = "biological_process",
                             GOGroupByDepthType groupingType =
                             GO_GROUP_BY_MINIMUM_DEPTH);

  // for each cutoff, group the highest functions which annotate that many or fewer genes such that there are no parent-child pairs in the group
  void groupFunctionsByCutoff(map< string, map< unsigned int, set< GOFunction * > > > &functionsByCutoff, const vector< unsigned int > &cutoffs, const map< BioFunction, unsigned int > &geneCounts, bool allowDescendants = false);
  // group functions using their minimum depths such that there are no parent-child pairs in a group
  // this version is needed because the original doesn't prevent parent-child pairs and only considers one category at a time
  void groupFunctionsByDepth2(map< string, map< unsigned int, set< GOFunction * > > > &functionsByDepth, bool allowDescendants = false);

  // group functions using their minimum depths such that there are no
  // parent-child pairs in a group this version is needed because i
  // want to return BioFunctions rather than pointers to GOFunctions.
  void groupFunctionsByDepth2(map< string, map< unsigned int, set< BioFunction > > > &functionsByDepth, bool allowDescendants = false);
  ~GeneOntology();

private:

  bool _boolRelationship;
  string _strRelationshipFile;
  
  void _init()
    {
      _functionCategoryAbbrevs["p"] = "biological_process";
      _functionCategoryAbbrevs["c"] = "cellular_component";
      _functionCategoryAbbrevs["f"] = "molecular_function";
      _legalFunctionCategories.insert("biological_process");
      _legalFunctionCategories.insert("cellular_component");
      _legalFunctionCategories.insert("molecular_function");
    }
  
  void addIdx(unsigned int id,GOFunction* obj);
  void buildIndices();

  bool _isLegalFunctionCategory(string name, string &expandedName) const
    {
      expandedName = name;
      if (_legalFunctionCategories.end() !=
          _legalFunctionCategories.find(name))
        // name is in _legalFunctionCategories
        return(true);
      // name may be an abbreviation for a legal category
      map< string, string >::const_iterator aitr =
        _functionCategoryAbbrevs.find(name);
      if (_functionCategoryAbbrevs.end()
          != (aitr = _functionCategoryAbbrevs.find(name)))
        // name is a legal abbreviation.
        {
          expandedName = aitr->second;
          return(true);
        }
      
      
      return(false);
    }


private:

  string default_namespace;
  vector<GOFunction*> functions;
  
  //represent the integerId->function mapping as a vector since
  //it is bounded in size generally linear in count, and
  //vector is instant for lookups which are going to
  //happen alot with Ids. we can use a map but a vector
  //is faster at the cost of using extra space.
  vector<GOFunction*> idxId;
  map<string,GOFunction*> idxName;

  // various constants associated with GO.
  map< string, string > _functionCategoryAbbrevs;
  set< string > _legalFunctionCategories;
  set< string > _knownFunctionIds;
  
  bool _topologicalSortComputed;
  vector< GOFunction* > _topologicallySortedFunctions;
  
  // the first key (a string) is the category. the second
  // key (unsigned int) is either the number of annotations of a
  // function or its minimum depth.
  map< string, map< unsigned int, set< GOFunction * > > > _functionsByCutoff;
  map< string, map< unsigned int, set< GOFunction * > > > _functionsByDepth;
  
                
};

template < typename InputIterator, typename OutputIterator, typename Extractor >
void GeneOntology::computeMostSpecificFunctionsTemplated(
  InputIterator functionsBegin, InputIterator functionsEnd,
  OutputIterator specificFunctionsBegin,
  Extractor _extractor)
{
  bool isAncestor;
  InputIterator fitr1, fitr2;
  GOFunction *func1, *func2;
  // compare all pairs of functions to see if one is an ancestor of the other.
  for (fitr1 = functionsBegin; fitr1 != functionsEnd; fitr1++)
    {
      // assume fitr1 is not an ancestor of any other function.
      isAncestor = false;
      func1 = getFunctionById(_extractor(*fitr1));
      fitr2 = fitr1;
      fitr2++;
      for (; fitr2 != functionsEnd; fitr2++)
        {
          func2 = getFunctionById(_extractor(*fitr2));
          if (func2->isAncestor(func1))
          // *fitr1 is an ancestor of *fitr2. i do not need to make
          // any more checks for *fitr1.
          {
            isAncestor = true;
            break;
          }
        }
      if (!isAncestor)
        // fitr1 is not an ancestor of any other functions. add it to
        // specificFunctionsBegin.
        {
          *specificFunctionsBegin = *fitr1;
          ++specificFunctionsBegin;
        }
    }
}

#endif // _GO_H