old-annotations.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 _ANNOTATIONS_H
#define _ANNOTATIONS_H

#include <iostream>
#include <fstream>
#include <iterator>
#include <map>
#include <set>
#include <string>

using namespace std;

#include "boost/multi_index_container.hpp"
#include "boost/multi_index/member.hpp"
#include "boost/multi_index/mem_fun.hpp"
#include "boost/multi_index/ordered_index.hpp"
#include <boost/multi_index/composite_key.hpp>

using boost::multi_index_container;
using namespace boost::multi_index;

// for is_any_of.
#include "boost/algorithm/string/classification.hpp"
// for string split.
#include "boost/algorithm/string.hpp"

#include "enrichment.h"
#include "myiterator.h"

#include "biofunction.h"
#include "GO.h"

// also declared in graph.h
typedef double MyNT;

//
// there is no good reason old-annotations.h should have this
// declaration. but MyAnnotations::computeEnrichmentsGenGO() has an
// argument of type GenGOParametersOptimisationType. i cannot #include
// gengo.h in this file since gengo.h refers to MyAnnotations.
//
enum GenGOParametersOptimisationType { GENGO_PARAMETER_OPTIMISATION_NONE = 0,
                                       GENGO_PARAMETER_OPTIMISATION_LOOP,
                                       GENGO_PARAMETER_OPTIMISATION_LEARN,
                                       // loop over many values and learn for each value.
                                       GENGO_PARAMETER_OPTIMISATION_LOOP_AND_LEARN,
};




enum MyGainAnnotationType { NOT_ANNOTATED_STATE = -1, HYPOTHETICAL_STATE,
                        ANNOTATED_STATE};

inline MyGainAnnotationType makeAnnotationType(string c)
{
  if ("-1" == c)
    return(NOT_ANNOTATED_STATE);
  if ("0" == c)
    return(HYPOTHETICAL_STATE);
  if ("1" == c)
    return(ANNOTATED_STATE);
  cerr << "Unknown annotation state" << c << endl;
  exit(-1);
}

extern map< string, string > GO_UNKNOWN_FUNCTIONS;

//     *   TAS / IDA
//     *
//           o IMP = / IGI / IPI
//           o ISS / IEP
//           o NAS
//           o
//                 + = IEA

struct GOEvidenceCodes
{
private:
  map< string, MyNT >_weights;
public:

  GOEvidenceCodes()
    {
        /*
      _weights["TAS"] = 0.9;
      _weights["IDA"] = 0.9;
      _weights["IMP"] = 0.7;
      _weights["IGI"] = 0.7;
      _weights["IPI"] = 0.7;
      _weights["IEP"] = 0.6;
      _weights["ISS"] = 0.6;
      _weights["IC"] = 0.5;
      _weights["NAS"] = 0.4;
      _weights["RCA"] = 0.4;
      _weights["IEA"] = 0.1;
      */
      
      // Brett Tyler's recommended weightings:
      //
      // "Note that TAS is weak because the nature of the evidence is not
      // determined. GO considers ISS and IEP to be special cases of RCA so
      // all should have similar weightings. IC is used very conservatively so
      // is usually reliable."
      _weights["IDA"] = 0.9;
      _weights["IMP"] = 0.9;
      _weights["IGI"] = 0.9;
      _weights["IC"] = 0.8;
      _weights["IPI"] = 0.6;
      _weights["IEP"] = 0.6;
      _weights["ISS"] = 0.6;
      _weights["RCA"] = 0.6;
      _weights["TAS"] = 0.6;
      _weights["NAS"] = 0.3;
      _weights["IEA"] = 0.2;
    }
  
  MyNT getWeight(string code) const
    {
      static set< string > _unknownCodes;
      map< string, MyNT >::const_iterator witr = _weights.find(code);
      if (_weights.end() == witr)
        {
          if (_unknownCodes.end() == _unknownCodes.find(code))
            {
              cerr << "GOEvidenceCodes::getWeight(): unknown evidence code "
                   << code << ", returning weight of 0" << endl;
              _unknownCodes.insert(code);
            }
          return(0);
        }
      
      return(witr->second);
    }

  bool hasWeight(string code) const
    {
      map< string, MyNT >::const_iterator witr = _weights.find(code);
      return(_weights.end() != witr);
    }
  

  bool isExperimental(string code) const
    {
      if ("TAS" == code)
        return(true);
      if ("IDA" == code)
        return(true);
      if ("IMP" == code)
        return(true);
      if ("IGI" == code)
        return(true);
      if ("IPI" == code)
        return(true);
      if ("IEP" == code)
        return(true);
      if ("ISS" == code)
        return(false);
      if ("IC" == code)
        return(true);
      if ("NAS" == code)
        return(true);
      if ("RCA" == code)
        return(false);
      if ("IEA" == code)
        return(false);
    }

  void checkEvidenceCodeWeights(const map< string, unsigned int > &ecUsageCounts, ostream &logStream);


  void readEvidenceCodeWeights(string ecwFile);
};

// read functions from file. each line contains one function.
template < class container >
void readFunctions(string file, container &functions)
{
  cout << "Reading items from " << file << endl;
  ifstream ifstr(file.c_str());
  if (!ifstr)
    {
      cerr << "Cannot open file " << file << endl;
      exit(-1);
    }
  string func, line;
  vector< string > items;
  while (!ifstr.eof())
    {
      getline(ifstr, line);
      // ignore comments.
      if ('#' == line[0])
        continue;
      items.clear();
      boost::split(items, line, boost::is_any_of("\t"));
      if (items.size() == 0)
        continue;
      func = items[0];
      if (func=="")
        continue;
      functions.insert(func);
    }
  ifstr.close();
  cout << "Read " << functions.size() << " items." << endl;
}
 
// tell doxygen to skip
struct annotation_gene_function_evidencecode_index {};
struct annotation_gene_function_index {};
struct annotation_gene_functiontype_index {};
struct annotation_gene_index {};
struct annotation_function_index {};
//struct algorithm_index {};

struct annotation_function_gene_index {};






struct FunctionalAnnotation
{
public:
  // these need to be public for SetOfAnnotations.
  string _gene;
  BioFunction _function;
  string _evidenceCode;
  MyGainAnnotationType _annotationStatus;
  // must be false when _annotationStatus is not ANNOTATED_STATE.
  bool _mostSpecificAnnotation;
  
public:
  FunctionalAnnotation(string g, string f, string c, string e,
                       MyGainAnnotationType status, bool msa = false)
      : _gene(g), _function(c, f), _evidenceCode(e),
        _annotationStatus(status),_mostSpecificAnnotation(msa)
    {}
  FunctionalAnnotation(string g, const BioFunction &f, string e,
                       MyGainAnnotationType status, bool msa = false)
      : _gene(g), _function(f), _evidenceCode(e), _annotationStatus(status),
        _mostSpecificAnnotation(msa)
    {}

  string getFunctionType() const
    {
      return(_function.getCategory());
    }
  
};


// // extracts a FunctionalAnnotation's function type. See "User-defined
// // key extractors" in
// // http://www.boost.org/doc/libs/1_36_0/libs/multi_index/doc/tutorial/key_extraction.html
// struct FunctionalAnnotationFunctionType
// {
//   // result_type typedef required by Key Extractor concept
//   typedef string result_type; 
  
//   string operator()(const FunctionalAnnotationFunctionType& fa) const // operator() must be const
//   {
//     return fa._function.getCategory();
//   }
// };


typedef multi_index_container<
  // the class whose instances we will store in the multi_index_container.
  FunctionalAnnotation,
  // specify all the indices now.
  indexed_by<

    // Define multi-column primary key.  Order first by
    // FunctionalAnnotation::_gene, then by
    // FunctionalAnnotation::_function, and finally by
    // FunctionalAnnotation::_evidenceCode. Use
    // annotation_gene_function_evidencecode_index as the tag for this
    // index. This index is unique since a specific (function, gene)
    // pair may occur exactly once. Hence, it serves as the primary key for the multi_index.
    ordered_unique<
      // name by which to refer to this index.
      tag< annotation_gene_function_evidencecode_index >,
      // specify the composite key now.
      composite_key<
        // the class whose members form the composite key. typically, it is the same class as the one we are storing in the multi_index.
        FunctionalAnnotation,
        // specify how to extract each element of the key. the first
        // item says that we want a member of FunctionalAnnotation,
        // its type is string, and it can be accessed through the
        // _gene member variable.
        member< FunctionalAnnotation, string, &FunctionalAnnotation::_gene>,
        member< FunctionalAnnotation, BioFunction, &FunctionalAnnotation::_function >,
        member< FunctionalAnnotation, string, &FunctionalAnnotation::_evidenceCode >
        >
//       composite_key_compare<
//       std::less<std::string>, // genes sorted by default
//       std::less<std::string>, // functions sorted by default
//       std::less<std::string>  // evidencecodes sorted by default
//       >
    >,
    
    // specify other non-primary indices now.
    //
    // there can be many evidence codes for a gene-function pair.
    ordered_non_unique<
    tag< annotation_gene_function_index >,
    composite_key<
    FunctionalAnnotation,
    member< FunctionalAnnotation, string, &FunctionalAnnotation::_gene >,
    member< FunctionalAnnotation, BioFunction, &FunctionalAnnotation::_function >
  >,
    composite_key_compare<
      std::less<std::string>, // genes sorted by default
      std::less< BioFunction >  // functions sorted by default
  >
  >,

    // Order first by FunctionalAnnotation::_function and then by
    // FunctionalAnnotation::_gene. Use annotation_function_gene_index
    // as the tag for this index. This index is non-unique since a
    // specific (function, gene) pair may occur multiple times.
    ordered_non_unique<
    tag< annotation_function_gene_index >,
    composite_key<
    FunctionalAnnotation,
    member< FunctionalAnnotation, BioFunction, &FunctionalAnnotation::_function >,
    member< FunctionalAnnotation, string, &FunctionalAnnotation::_gene >
      >,
      composite_key_compare<
        std::less< BioFunction >,  // functions sorted by default
        std::less<std::string> // genes sorted by default
        >
      >,

    
    // Order first by FunctionalAnnotation::_gene and then the value
    // returned by the member function
    // FunctionalAnnotation::getFunctionType(). Use
    // annotation_gene_functiontype_index as the tag for this
    // index. This index is non-unique since a specific (gene,
    // functiontype) pair may occur multiple times.
    ordered_non_unique<
      tag< annotation_gene_functiontype_index >,
      composite_key<
        FunctionalAnnotation,
        member< FunctionalAnnotation, string, &FunctionalAnnotation::_gene >,
        // see http://www.boost.org/doc/libs/1_36_0/libs/multi_index/doc/tutorial/key_extraction.html#const_mem_fun
        const_mem_fun< FunctionalAnnotation, string, &FunctionalAnnotation::getFunctionType >
        >
      >,

    // sort by less<string> on _gene
    ordered_non_unique< tag< annotation_gene_index >,
                        member< FunctionalAnnotation, std::string,
                                &FunctionalAnnotation::_gene > >,
    
    // sort on _function
    ordered_non_unique< tag< annotation_function_index >,
                        member< FunctionalAnnotation, BioFunction,
                                &FunctionalAnnotation::_function > >//,
    
    > // matches indexed_by 
  > SetOfFunctionalAnnotations;

class MyAnnotations
{

private:
//  typedef map< string, map< string, map< string, MyGainAnnotationType > > >
  typedef set< string > BioFunctionSet;
  typedef map< string, BioFunctionSet > MyGenesToBioFunctionTypes;
  typedef map< string,  MyGenesToBioFunctionTypes > MyGenesToFunctions;
  // first "string" is the function type, second "string" is the
  // function, third "string" is the gene/protein.
//  typedef map< string, map< string, MyGainAnnotationType > > BioFunctionToGenes;
  typedef set< string > MyGeneSet;
  typedef map< string, MyGeneSet > BioFunctionToGenes;
  typedef map< string, BioFunctionToGenes > BioFunctionTypesToGenes;

  typedef BioFunctionToGenes::iterator _BioFunctionsIterator;
  typedef BioFunctionToGenes::const_iterator _BioFunctionsConstIterator;
  typedef BioFunctionTypesToGenes::iterator _NewBioFunctionsIterator;

  // multi_index related typedefs to make accessing various indices easier.
  typedef SetOfFunctionalAnnotations::index< annotation_function_index >::type FunctionIndex;
  typedef SetOfFunctionalAnnotations::index< annotation_function_gene_index >::type FunctionGeneIndex;
  typedef SetOfFunctionalAnnotations::index< annotation_gene_index >::type GeneIndex;

  
  SetOfFunctionalAnnotations _allAnnotations;
  
  SetOfFunctionalAnnotations _allNonAnnotations;
  
  // default constructor will initialise weights, which can be overridden by MyAnnotations::readEvidenceCodeWeights().
  GOEvidenceCodes _evidenceCodes;
  // the number of times each evidence code is used in an annotation.
  map< string, unsigned int > _evidenceCodeCountsInAllAnnotations;
  
  // list of genes annotated by at least one function of a given type.
//  map< string, map< string, MyGainAnnotationType > > functionTypeAnnotatedGenes;
  map< string, set< string > > functionTypeAnnotatedGenes;

  // list of genes not annotated by at any function of a given type.
//  map< string, map< string, MyGainAnnotationType > > functionTypeNonAnnotatedGenes;
  map< string, set< string > > functionTypeNonAnnotatedGenes;
  
  // for each function, list of genes annotating it.
  BioFunctionTypesToGenes functionsAnnotatedGenes;

  // for each function, list of genes non-annotated/hypothetical wrt function.
  BioFunctionTypesToGenes functionsNonAnnotatedGenes;

  
MyGenesToFunctions geneAnnotatedFunctions;
  MyGenesToFunctions geneNonAnnotatedFunctions;

  // for each function, the number of genes annotated with the function. 
  map< string, unsigned int > numFunctionAnnotatedGenes;
  // for each function, the number of genes
  // non-annotated/hypothetical with the function.
  map< string, unsigned int > numFunctionNonAnnotatedGenes;
  
  map< string, MyNT > functionFrequencies;

  // GO related variables, has true path rule been applied.
  GeneOntology *goDAG;
  bool _truePathRuleAppliedUpward;
  bool _truePathRuleAppliedDownward;
  bool _truePathRuleAppliedSideways;

  map< string, set< string > > _functionsWithOverlappingAnnoations;
  
  bool _overlappingFunctionsFlag;
  
public:  

  class NewBioFunctionsIterator;
  friend class NewBioFunctionsIterator;
  
  class BioFunctionsIterator: public  MyIterator< string >
  {
  private:
    _BioFunctionsIterator start;
    _BioFunctionsIterator end;
    _BioFunctionsIterator current;
    
  public:
    BioFunctionsIterator(const _BioFunctionsIterator _start,
                         const _BioFunctionsIterator _end)
        : start(_start), end(_end), current(_start)
      {}
    
    virtual bool hasNext()
      {
    return(current != end);
      }
    
    string &next()
      {
        // BUBUGBUG UGLY HACK. I have defined MyIterator::next() to return
        // a T& so that when T is a class, i can return a ref. but when T
        // is a string, i get compilation errors without this const_cast.
        string &_next(const_cast< string & >((*current).first));
        current++;
        return(_next);
      }
  };

  class BioFunctionsConstIterator: public  MyIterator< string >
  {
  private:
    _BioFunctionsConstIterator start;
    _BioFunctionsConstIterator end;
    _BioFunctionsConstIterator current;
    
  public:
    BioFunctionsConstIterator(_BioFunctionsConstIterator _start,
                              _BioFunctionsConstIterator _end)
        : start(_start), end(_end), current(_start)
      {}
    
    virtual bool hasNext()
      {
        return(current != end);
      }
    
    string &next()
      {
        // BUBUGBUG UGLY HACK. I have defined MyIterator::next() to return
        // a T& so that when T is a class, i can return a ref. but when T
        // is a string, i get compilation errors without this const_cast.
        string &_next(const_cast< string & >((*current).first));
        // string & _next((*current).first);
        current++;
        return(_next);
      }
  };
  
class BioFunctionsIteratorMultiIndex: public  MyIterator< BioFunction >
{
private:
  // the reference to the set of all annotations passed in the constructor.
  const SetOfFunctionalAnnotations & _allAnnotations;
  // the reference to the index over the functions in _allAnnotations.
  const FunctionIndex &_functionIndex;
  FunctionIndex::iterator _startItr;
  FunctionIndex::iterator _endItr;
  // the pointer to the current function.
  FunctionIndex::iterator _currentItr;
  // store the id of the current function.
  BioFunction _currentFunction;

  
public:
  BioFunctionsIteratorMultiIndex(const SetOfFunctionalAnnotations &_anns)
      : _allAnnotations(_anns), _functionIndex(_allAnnotations.get< annotation_function_index >()),
        _startItr(_allAnnotations.get< annotation_function_index >().begin()),
        _endItr(_allAnnotations.get< annotation_function_index >().end()), _currentItr(_startItr), 
        _currentFunction()
    {
      // sometime _startItr might be equal to _endItr, so I am being careful about initialising _currentFunction.
      if (_startItr != _endItr)
        _currentFunction = _currentItr->_function;
    }
  virtual bool hasNext()
  {
    return(_currentItr != _endItr);
  }
  BioFunction &next()
  {
    _currentFunction = _currentItr->_function;
    _moveCurrentItr();
    return(_currentFunction);
  }
private:      
  // move _currentItr to point to the start of the next function.
  void _moveCurrentItr()
    {
      _currentItr = _functionIndex.upper_bound(_currentFunction);
    }
  
};
  
  
public:
  MyAnnotations()
    : functionTypeAnnotatedGenes(), functionTypeNonAnnotatedGenes(),
         functionsAnnotatedGenes(), functionsNonAnnotatedGenes(),
         geneAnnotatedFunctions(),
         geneNonAnnotatedFunctions(),
         numFunctionAnnotatedGenes(),
      numFunctionNonAnnotatedGenes(),
      _evidenceCodes(), _evidenceCodeCountsInAllAnnotations(),
      functionFrequencies(),
      goDAG(NULL), _truePathRuleAppliedUpward(false), _truePathRuleAppliedDownward(false),
      _truePathRuleAppliedSideways(false), _overlappingFunctionsFlag(false)
  {
    // set up static variables.
    _initialise();
  }

  virtual ~MyAnnotations()
  {}

  void setOverlappingFunctionsFlag()
  {
     _overlappingFunctionsFlag = true;
  }

  void addAnnotation(string gene, const BioFunction &function,
                     string evidenceCode,
                     MyGainAnnotationType status = ANNOTATED_STATE,
                     bool msa = false);

  void addAnnotation(string gene, string function, string funcType,
                     string evidenceCode,
                     MyGainAnnotationType status = ANNOTATED_STATE,
                     bool msa = false);
  

  

  void applyTruePathRule(GeneOntology &go, bool applyUpward = true, bool applyDownward = true, bool applySideways = true);


  bool checkTruePathRule(GeneOntology &go, ostream *ostr = NULL,
                         bool *annotationsOK = NULL, bool *unknownOK = NULL);// const;


  void computeDifference(MyAnnotations &other, MyAnnotations &difference,
                         bool restrictToCommonGenes = false);

  void computeEnrichments(const set< string > &geneSet,
                          const GeneOntology *go,
                          vector<EnrichmentRecord< string, string> > &rejalts);
  
  void computeEnrichmentsGenGO(const set< string > &geneSet,
                               vector<EnrichmentRecord< string, string> > &rejalts,
                               const GeneOntology *go = NULL,
                               GenGOParametersOptimisationType optType = GENGO_PARAMETER_OPTIMISATION_NONE
                               ) const;

  void computeEnrichmentsGenGO(const set< string > &geneSet,
                               string functionType,
                               vector<EnrichmentRecord< string, string> > &rejalts,
                               const GeneOntology *go = NULL,
                               GenGOParametersOptimisationType optType = GENGO_PARAMETER_OPTIMISATION_NONE
                               ) const;

  MyNT computeFunctionFrequency(string funcType, string func);

  // compute various statistics.
  void computeStatistics();
  
  // count the number of genes each function annotates (NOT the number of annotations)
  void computeGeneCounts(map< BioFunction, unsigned int > &geneCounts);
  
  void copyAnnotationsForFunctions(const set< BioFunction > &functions, MyAnnotations &annotations);
  

  virtual MyAnnotations::BioFunctionsIterator functions(string type)
    {
      MyAnnotations::BioFunctionsIterator itr(functionsAnnotatedGenes[type].begin(),
                                              functionsAnnotatedGenes[type].end());
      return(itr);
    }

  virtual MyAnnotations::BioFunctionsIteratorMultiIndex functionsMultiIndex() const
    {
      // MyAnnotations::BioFunctionsIteratorMultiIndex itr(_allAnnotations.begin(), _allAnnotations.end(), _allAnnotations);
      MyAnnotations::BioFunctionsIteratorMultiIndex itr(_allAnnotations);
      return(itr);
    }
  
  set< string > functionTypes() const
    {
      set< string > types;
      BioFunctionTypesToGenes::const_iterator itr;
      for (itr = functionsAnnotatedGenes.begin(); itr != functionsAnnotatedGenes.end(); itr++)
        types.insert((*itr).first);
      return(types);
    }

  MyNT getAnnotatedFunctionWeight(string gene, string func, string type);

  MyNT getAnnotatedFunctionWeightProbabilisticOr(string gene, string func, string type);
  
  map< string, MyNT > getAnnotatedFunctionWeights(string gene, string func, string type);

  // note: this method is used in Reporter::printECWeightedCurveDataFromCV().
  map< string, MyNT > getAnnotatedWeights(string gene, string type);

  void getAnnotatedGenes(string functype, set< string >& annotatedGenes);

  void getAnnotations(string type, map< string, set< string > > &annotations) const;
  
  void getAnnotationsByGene(string type, map< string, set< string > > &annotations) const;
  
  virtual void getAnnotations(string type, string function,
                              set< string >  &annotations) const;
  
  void getAnnotations(const GOFunction * function, set< string > &annotations) const;

  void getAnnotations(const set< GOFunction *>  &functions,
                      set< string >  &annotations) const;

  virtual void getAnnotationsForGene(string type, string gene,
                                     set< string >  &annotations) const;


  void getEvidenceCodeCounts(map< string, unsigned int > &ecCounts) const
    {
      ecCounts = _evidenceCodeCountsInAllAnnotations;
    }
  
  void getEvidenceCodeCounts(string funcType, string funcId, map< string, unsigned int > &ecCounts) const;
  
  
  
  MyNT getFunctionFrequency(string function) const
    {
      map< string, MyNT >::const_iterator itr;
      itr = functionFrequencies.find(function);
      if (itr == functionFrequencies.end())
        return(0);
      return((*itr).second);
    }

  string getFunctionType(string functionName) const;
  

  bool hasExperimentallyAnnotatedFunction(string gene, string category) const
    {
      // TODO is there a better way to do this perhaps by using a gene-evidencecode index?
      const GeneIndex &gi = _allAnnotations.get< annotation_gene_index >();
      GeneIndex::const_iterator start = gi.lower_bound(gene), end = gi.upper_bound(gene);
      for (GeneIndex::const_iterator gitr = start; gitr != end; gitr++)
        {
          if (_evidenceCodes.isExperimental(gitr->_evidenceCode) && gitr->_function.getCategory() == category)
            return (true);
        }
      return (false);
    }
  
        
  bool hasAnnotatedFunction(string gene) const
    {
      const GeneIndex &gi = _allAnnotations.get< annotation_gene_index >();
      return (gi.end() != gi.find(gene));
      //return(geneAnnotatedFunctions.end() !=
      //       geneAnnotatedFunctions.find(gene));
    }

  bool hasAnnotatedFunction(string gene, string category) const
    {
      return(hasAnnotatedFunction(gene) &&
             // how i can check in this category.
             (geneAnnotatedFunctions.find(gene)->second.end() !=
              geneAnnotatedFunctions.find(gene)->second.find(category)));
    }
        
  // method's name used to be hasNoFunction().
   bool isNotKnownToBeAnnotatedWithFunction(string gene, string func, string type)
    {
      if (!_truePathRuleAppliedDownward)
        // use the explicit check.
        return(isNotKnownToBeAnnotatedWithFunctionSansTruePathRule(gene, func, type, *goDAG));
    
      if(_overlappingFunctionsFlag)
      {
          if (haveOverlappingFunction(gene, func, type))
             return true;      
      }


      // since _truePathRuleAppliedDownward is false, use
      // geneNonAnnotatedFunctions to answer the query.
      MyGenesToFunctions::const_iterator gfItr;
      MyGenesToBioFunctionTypes::const_iterator gbftItr;
      // gene may have no function stored at all!
      if (geneAnnotatedFunctions.end() ==
          (gfItr = geneAnnotatedFunctions.find(gene)))
        return(true);
      // gene may have no function stored for type.
      if (gfItr->second.end() == (gbftItr = gfItr->second.find(type)))
        return(true);
      // function stored for type is unknown function.
      if (gbftItr->second.end() !=
          gbftItr->second.find(GO_UNKNOWN_FUNCTIONS[type]))
        return(true);
      // gene is non-annotated with respect to function, i.e. an
      // ancestor of function is the most specific function
      // annotating gene.
      if ((geneNonAnnotatedFunctions[gene][type].end() !=
           geneNonAnnotatedFunctions[gene][type].find(func))
//           &&
//           // this check is redundant because i put a func in
//           // geneNonAnnotatedFunctions only if the state is
//           // HYPOTHETICAL_STATE.
//           (HYPOTHETICAL_STATE ==
//            geneNonAnnotatedFunctions[gene][type][func])
          )
        return(true);

      // gene may have no function stored at all!
      if (geneAnnotatedFunctions.end() ==
          geneAnnotatedFunctions.find(gene))
        return(true);
      // gene may have no function stored for type.
      if (geneAnnotatedFunctions[gene].end() ==
         geneAnnotatedFunctions[gene].find(type))
        return(true);
      // function stored for type is unknown function.
      if (geneAnnotatedFunctions[gene][type].end() !=
          geneAnnotatedFunctions[gene][type].find(GO_UNKNOWN_FUNCTIONS[type]))
        return(true);
      // gene is non-annotated with respect to function, i.e. an
      // ancestor of function is the most specific function
      // annotating gene.
      if ((geneNonAnnotatedFunctions[gene][type].end() !=
           geneNonAnnotatedFunctions[gene][type].find(func))
//           &&
//           // this check is redundant because i put a func in
//           // geneNonAnnotatedFunctions only if the state is
//           // HYPOTHETICAL_STATE.
//           (HYPOTHETICAL_STATE ==
//            geneNonAnnotatedFunctions[gene][type][func])
          )
        return(true);


      return(
        // gene may have no function stored at all!
        (geneAnnotatedFunctions.end() ==
         geneAnnotatedFunctions.find(gene)) || 
        // gene may have no function stored for type.
        (geneAnnotatedFunctions[gene].end() ==
         geneAnnotatedFunctions[gene].find(type)) ||
        // function stored for type is unknown function.
        (geneAnnotatedFunctions[gene][type].end() !=
         geneAnnotatedFunctions[gene][type].find(GO_UNKNOWN_FUNCTIONS[type])) ||
        // gene is non-annotated with respect to function, i.e. an
        // ancestor of function is the most specific function
        // annotating gene.
        ((geneNonAnnotatedFunctions[gene][type].end() !=
          geneNonAnnotatedFunctions[gene][type].find(func))
         //          &&
         //          // this check is redundant because i put a func in
         //          // geneNonAnnotatedFunctions only if the state is
         //          // HYPOTHETICAL_STATE.
         //          (HYPOTHETICAL_STATE ==
         //           geneNonAnnotatedFunctions[gene][type][func])
         )
        );
    }

  bool isNotKnownToBeAnnotatedWithFunctionSansTruePathRule(string gene, string func, string type,
                                                         GeneOntology &goDAG);


  bool hasAnnotatedGene(string func, string type) const
    {
      BioFunctionTypesToGenes::const_iterator itr = functionsAnnotatedGenes.find(type);
      if (functionsAnnotatedGenes.end() == itr)
        // this type has no function.
        return(false);
      BioFunctionToGenes::const_iterator gitr = itr->second.find(func);
      if (itr->second.end() == gitr)
        // the function does not annotate any gene.
        return(false);
      // check if there indeed are any genes stored at gitr.
      if (gitr->second.empty())
        return(false);
      return(true);
    }
  
  bool hasAnnotatedGene(const BioFunction &function) const
    {
      return(hasAnnotatedGene(function.getId(), function.getCategory()));
    }
  
  
  // return true if gene1 and gene2 have the given function.
  bool haveSameFunction(string gene1, string gene2, const BioFunction &function);
  // return true if gene1 and gene2 are both annotated with at least
  // one of the functions in the vector functionVector.
  bool haveSameFunction(string gene1, string gene2,
                                       const vector< BioFunction >& functionVector);
  template < class InputIterator >
  bool haveSameFunction(string gene1, string gene2,
                                       InputIterator first, InputIterator last);
  // return true if gene1 and gene2 have the same function of the given type.
  bool haveSameFunction(string gene1, string gene2, string type);
  // return true if gene1 and gene2 have the same function (of any type).
  bool haveSameFunction(string gene1, string gene2);

  // return true if the annotation of gene1 overlapps with function.
  bool haveOverlappingFunction(string gene1, string function, string type);


  
  //
  // WARNING: this method works correctly only with the multi_index
  // structure, since the multi_index stores whether an annotation is
  // the most specific annotation or not.
  bool isAnnotatedWithFunction(string gene, string func, string type, bool beforeTruePathRule = false);
  
 
  bool isAnnotatedWithFunctionExperimentally(string gene, string func, string type);

  bool isUnknownFunction(string func, string type) const
    {
      map< string, string >::const_iterator itr = GO_UNKNOWN_FUNCTIONS.find(type);
      return((GO_UNKNOWN_FUNCTIONS.end() != itr)
             && (func == itr->second));
    }

  //void keepAnnotationsForFunctions(set< string > &functionSet);

  void keepAnnotationsForGenes(set< string > &geneSet);

  
  unsigned int numAnnotations(string category = "");
  
  // use to be numNonAnnotations.
  unsigned int numUnknownAnnotations(string category = "");
  
  unsigned int numAnnotatedGenes(string funcType, string function)
    {
      return(functionsAnnotatedGenes[funcType][function].size());
//    return(_allAnnotations.get< annotation_function_index >().count(function));
    }

  unsigned int numAnnotatedGenes(string funcType = "")
  {
    if ("" != funcType)
      return(functionTypeAnnotatedGenes[funcType].size());
    set< string > funcTypes = functionTypes();
    set< string >::iterator ftItr;
    unsigned int sum = 0;
    for (ftItr = funcTypes.begin(); ftItr != funcTypes.end(); ftItr++)
      sum += functionTypeAnnotatedGenes[*ftItr].size();
      return(sum);
  }
  
  unsigned int numAnnotatingFunctions(string gene, string funcType)
  {
    return(geneAnnotatedFunctions[gene][funcType].size());
  }

  unsigned int numAnnotatingFunctions(string funcType = "")
    {
      if ("" != funcType)
        return(functionsAnnotatedGenes[funcType].size());
      set< string > funcTypes = functionTypes();
      set< string >::iterator ftItr;
      unsigned int sum = 0;
      for (ftItr = funcTypes.begin(); ftItr != funcTypes.end(); ftItr++)
        sum += numAnnotatingFunctions(*ftItr);
      return(sum);
    }
  
  unsigned int numNonAnnotatedGenes(string funcType)
  {
#ifdef DEBUG
    cout << "function type is " << funcType
         << ", unknown function is " << GO_UNKNOWN_FUNCTIONS[funcType]
         << ", #non-annotated genes is "
         << functionsNonAnnotatedGenes[funcType][GO_UNKNOWN_FUNCTIONS[funcType]].size()
         << " or " << functionTypeNonAnnotatedGenes[funcType].size()
         << endl;
    
#endif    
    return(functionTypeNonAnnotatedGenes[funcType].size());
  }

  unsigned int numTotalNonAnnotatedGenes(string function, string funcType)
  {
    return(numNonAnnotatedGenes(funcType) + 
           numNonAnnotatedGenes(funcType, function));
  }

  // return only how many genes are non-annotated wrt function.
  unsigned int numNonAnnotatedGenes(string funcType, string function)
  {
    return(functionsNonAnnotatedGenes[funcType][function].size());
  }
  
  // was numNonAnnotatingFunctions.
  unsigned int numUnknownAnnotatingFunctions(string gene, string funcType)
  {
    return(geneNonAnnotatedFunctions[gene][funcType].size());
  }


  // a gene is negative wrt function if it is annotated in funcType
  // (this number will/should not include genes that are not
  // annotated by any function in funcType) and does not annotate
  // function and is also not non-annotated wrt function.
  unsigned int numNegativeGenes(string funcType, string function)
  {
    return(numAnnotatedGenes(funcType) - 
           numAnnotatedGenes(funcType, function) -
           numNonAnnotatedGenes(funcType, function)
           );
  }

  void print(string fileName);

  void printWithName(string fileName);

  void printNumAnnotationsPerFunction(ostream &numFunStream);
  
  void read(string fileName, string dataType = "unweighted",
            const map< string, set< string > > *nodeAliases = NULL,
            bool mostSpecificAnnotations = false,
            const set< string > *evidenceCodesToIgnore = NULL);


  void read(string fileName, unsigned int keySize);
  

  void readGO(string fileName, unsigned int altSymbolCol=0);

  void readEvidenceCodeWeights(string ecwFile)
    {
      _evidenceCodes.readEvidenceCodeWeights(ecwFile);
      _evidenceCodes.checkEvidenceCodeWeights(_evidenceCodeCountsInAllAnnotations, cerr);
    }

  void subtract(MyAnnotations &other, MyAnnotations &difference);


  void storeGODAG(GeneOntology &go)
    {
      goDAG = &go;
    }

  void readOverlappingFunctions(string fileName);
  
private:

  // remove genes and functions with no annotations.
  void _cleanUp();
  
  unsigned int _copyAnnotatedGenes(const BioFunction &source, const BioFunction &target);

  unsigned int _copyAnnotatedGenesToNonAnnotated(const BioFunction &source, const BioFunction &target,
                                                 bool force = false,
                                                 MyAnnotations *original = NULL);
  

  
  unsigned int _copyNonAnnotatedGenes(const BioFunction &source, const BioFunction &target);

// helper function for MyAnnotations::isAnnotatedWithFunction() that
// uses the multi-level map.
bool _isAnnotatedWithFunctionMap(string gene, string func, string type);

// helper function for MyAnnotations::isAnnotatedWithFunction() that
// uses the multi-index. Read the documentation for MyAnnotations::isAnnotatedWithFunction to see the role of beforeTruePathRule.
bool _isAnnotatedWithFunctionMultiIndex(string gene, string func, string type, bool beforeTruePathRule = false);

  map< string, vector< string > > _functionsWithOverlappingAnnotations;

private:

  static void _initialise();
  void _compareFunctionTypes(string funcType1, string funcType2);
  bool _check();

};

  
  
  

#endif // _ANNOTATIONS_H