enrichment-algorithm.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/>.      *
 *                                                                        *
 **************************************************************************/

// Purpose: define classes that implement various algorithms for
// functional enrichment.

#ifndef _ENRICHMENT_ALGORITHM_H
#define _ENRICHMENT_ALGORITHM_H

#include "graph.h"
#include "old-annotations.h"

#include "../libexpression/point.h"
#include "constants.h"
// Test for GCC >= 4.3. from http://gcc.gnu.org/onlinedocs/cpp/Common-Predefined-Macros.html.
// See definition of GCC_VERSION macro in libutil/constants.h
#if GCC_VERSION >= 40300
// see http://gcc.gnu.org/gcc-4.3/changes.html
#include <tr1/unordered_set>
#include <tr1/unordered_map>
using namespace std::tr1;
#else
#include <ext/hash_map>
using namespace __gnu_cxx;
#define unordered_set hash_set;
#define unordered_map hash_map;
#endif

using namespace __gnu_cxx;

class LogHelper
{
private:
  mutable vector<double> _logFactCache;
public:
  LogHelper();
  double factorial(unsigned int) const;
  double factorialLog(unsigned int) const;
  double factorialEstimate(unsigned int) const;
  double factorialEstimateLog(unsigned int) const;
};


class FunctionalEnrichmentAlgorithm 
{
protected:
  set<MyNodeId> _clusterGenes;

  set<MyNodeId> _universeGenes;

  MyAnnotations& _annotations;

  const GeneOntology *_go;

  string _category;


  // This is declared as static because the LogHelper caches factorial logs
  // as they are computed. Now when any instance of FunctionalEnrichmentAlgorithm
  // is instantiated, it can share this cache (and the associated LogHelper).
  static LogHelper _logHelper;

  
public:
  FunctionalEnrichmentAlgorithm(const set< MyNodeId >& clusterGenes, const set< MyNodeId >& universeGenes,
                                MyAnnotations& annotations, const GeneOntology& go, string category);

  FunctionalEnrichmentAlgorithm(MyAnnotations& annotations, const GeneOntology& go, string category)
    : _annotations(annotations), _go(&go), _category(category) {}

  virtual ~FunctionalEnrichmentAlgorithm()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  virtual void computeEnrichmentForFunction(const string funcType, const string funcId) = 0;

  virtual void printEnrichmentResults(ostream &ostr) = 0;

};


class RandomFunctionNetworkBasedEnrichmentAlgorithm: public FunctionalEnrichmentAlgorithm
{
private:

  MyGraph _clusterGraph;
  MyGraph _universeGraph;

//  struct randomFunctionEnrichmentResult{
//    string functionId;
//    string functionType;
//    unsigned int v_f;
//    unsigned int c_f;
//    double lexComponentsByNodesPValue;
//    double lexComponentsByEdgesPValue;
//    double densityPValue;
//  };

  struct randomFunctionEnrichmentResult{
    string functionId;
    string functionType;
    unsigned int v_f;
    unsigned int c_f;
    double lexComponentsByNodesPValue;
    double lexComponentsByEdgesPValue;
    double analyticalHypPValue;
    double empiricalHypPValue;
    double densityPValue;
  };

  vector< randomFunctionEnrichmentResult > _enrichmentResults;

  map< string, vector< int > > _componentSizesByNodes;

  map< string, vector< int > > _componentSizesByEdges;

  map< string, double > _densities;

  map< string, pair< unsigned int, unsigned int > > _functionSizes;

  map< string, pair< int, int > > _lexComponentsByNodesRecords;

  map< string, pair< int, int > > _lexComponentsByEdgesRecords;

  map< string, pair< int, int > > _empiricalHypRecords;


  map< string, pair< int, int > > _densityRecords;

  unsigned int _numRandomGraphs;

  // a helper function for initializing _componentSizesByNodes, _densities
  void initialize(const set< string >& functionsToProcess);

public:
  RandomFunctionNetworkBasedEnrichmentAlgorithm(const MyGraph& clusterGraph, const MyGraph& universeGraph,
                 MyAnnotations& annotations, const GeneOntology& go, string category,
                 unsigned int numRandomGraphs);

  virtual ~RandomFunctionNetworkBasedEnrichmentAlgorithm()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);


  virtual void computeEnrichmentForFunction(const string funcType, const string funcId);

  void printEnrichmentResult(ostream &ostr, randomFunctionEnrichmentResult result);

  virtual void printEnrichmentResults(ostream &ostr);

};


class RandomUniverseNetworkBasedEnrichmentAlgorithm: public FunctionalEnrichmentAlgorithm
{
private:

  MyGraph _clusterGraph;
  MyGraph _universeGraph;
  MyGraph _randomGraph;

  struct randomUniverseEnrichmentResult{
    string functionId;
    string functionType;
    unsigned int v_f;
    unsigned int c_f;
    double lexComponentsByNodesPValue;
    double lexComponentsByEdgesPValue;
    double analyticalHypPValue;
    double empiricalHypPValue;
    double densityPValue;
  };

  vector< randomUniverseEnrichmentResult > _enrichmentResults;

  map< string, vector< int > > _componentSizesByNodes;

  map< string, vector< int > > _componentSizesByEdges;

  map< string, double > _densities;

  map< string, pair< unsigned int, unsigned int > > _functionSizes;

  map< string, pair< int, int > > _lexComponentsByNodesRecords;

  map< string, pair< int, int > > _lexComponentsByEdgesRecords;

  map< string, pair< int, int > > _empiricalHypRecords;

  map< string, pair< int, int > > _densityRecords;

  unsigned int _numRandomGraphs;

  unsigned int _numRandomizationIterations;

  void updateEnrichmentForFunction(string funcType, string funcId);

  // a helper function for initializing _componentSizesByNodes, _densities
  void initialize(const set< string >& functionsToProcess);

public:
  RandomUniverseNetworkBasedEnrichmentAlgorithm(const MyGraph& clusterGraph, const MyGraph& universeGraph,
                 MyAnnotations& annotations, const GeneOntology& go, string category,
                 unsigned int numRandomGraphs, unsigned int numRandomizationIterations);

  virtual ~RandomUniverseNetworkBasedEnrichmentAlgorithm()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  void computeEnrichmentForFunction(const string funcType, const string funcId);

  void printEnrichmentResult(ostream &ostr, randomUniverseEnrichmentResult result);

  virtual void printEnrichmentResults(ostream &ostr);

};


class HypergeometricAlgorithm: public FunctionalEnrichmentAlgorithm
{
public:
  struct hypergeometricEnrichmentResult{
    string functionId;
    string functionType;
    unsigned int v_f;
    unsigned int c_f;
    double hypergeometricPvalue;
  };

private:

  vector< hypergeometricEnrichmentResult > _hypergeometricEnrichmentResults;

public:

  HypergeometricAlgorithm(const set<string> &clusterGraph, const set<string> &universeSet,
                          MyAnnotations &annotations, const GeneOntology &go,
                          string category);


  HypergeometricAlgorithm(const MyGraph& clusterGraph, const MyGraph& universeGraph,
                          MyAnnotations& annotations, const GeneOntology& go,
                          string category);

  virtual ~HypergeometricAlgorithm()
    {}

  virtual void computeEnrichmentForFunction(const string funcType, const string funcId);

  virtual void printEnrichmentResults(ostream &ostr);

  static long double getHypergeometricPvalue(int v, int c, int v_f, int c_f);
};



class PAGEAlgorithm : public FunctionalEnrichmentAlgorithm
{

private:

  MyPointSet _universePoints;

  map< string, unsigned int > _indexMap;

  double _mean;
  double _stddev;

  double getValue(string gene)
    {
      return _universePoints[0][_indexMap[gene]];
    }

  struct enrichmentResult{
    string funcId;
    string funcType;
    unsigned int v_f;
    double zscore;
    double pval;
  };

  vector< enrichmentResult > _enrichmentResults;


public:

  PAGEAlgorithm(const MyPointSet &universePoints,
                MyAnnotations &annotations, const GeneOntology &go, string category);

  virtual ~PAGEAlgorithm()
    {}

  virtual void computeEnrichmentForFunction(const string funcType, const string funcId);

  virtual void printEnrichmentResults(ostream &ostr);
};



class GenGOAlgorithm: public FunctionalEnrichmentAlgorithm
{
private:

  double _p;
  double _q;
  double _penalty;

  struct FunctionInfo
  {
    set<MyNodeId> clusterGenes;
    int numNotClusterGenes;
  };
  map<string, FunctionInfo> _allFunctionInfo;
  vector<string> _activeFunctions;
  set<string> _notActiveFunctions;

  map<string, int> _clusterGenesToNumActiveFunctions;

  int _numClusterGenesWithSomeActiveFunction;

  int _numClusterGenesWithNoActiveFunction;

  int _numEdgesNotClusterGeneWithActiveFunction;

  int _numEdgesNotClusterGeneWithNotActiveFunction;

  void _activateFunction(string func);

  void _deactivateFunction(string func);

  double _calcLogLikelihood();

public:
  GenGOAlgorithm(const MyGraph& clusterGraph, const MyGraph& universeGraph,
                 MyAnnotations& annotations, const GeneOntology& go, string category,
                 double p, double q, double penalty);
  virtual ~GenGOAlgorithm()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  // not used
  virtual void computeEnrichmentForFunction(const string funcType, const string funcId)
    {}

  virtual void printEnrichmentResults(ostream &ostr);

};


class GroupHypAlgorithm: public FunctionalEnrichmentAlgorithm
{
private:
  MyGraph _universeGraph;
  MyGraph _clusterGraph;

  map< string, set<MyNodeId> > _functionToAnnotatedClusterGenes;
  map< string, set<MyNodeId> > _functionToAnnotatedNotClusterGenes;

  vector<string> _activeFunctions;
  set<string> _notActiveFunctions;

  map<string, int> _clusterGenesToNumActiveFunctions;
  map<string, int> _notClusterGenesToNumActiveFunctions;

  int _numClusterGenesWithSomeActiveFunction;
  int _numNotClusterGenesWithSomeActiveFunction;

  void _activateFunction(string func);

  void _deactivateFunction(string func);

  double _calcGroupHypPvalue() const;

public:
  GroupHypAlgorithm(const MyGraph& clusterGraph, const MyGraph& universeGraph,
                    MyAnnotations& annotations, const GeneOntology& go, string category);

  virtual ~GroupHypAlgorithm()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  // not used
  virtual void computeEnrichmentForFunction(const string funcType, const string funcId)
    {}

  virtual void printEnrichmentResults(ostream &ostr);

};




class EdgeGroupHyp: public FunctionalEnrichmentAlgorithm
{
private:
  MyGraph _universeGraph;
  MyGraph _clusterGraph;

  set<string> _clusterEdges;
  set<string> _notClusterEdges;


  map< string, set<string> > _functionToAnnotatedClusterGenes;

  map< string, set<string> > _functionToAnnotatedClusterEdges;
  map< string, set<string> > _functionToAnnotatedNotClusterEdges;

  vector<string> _activeFunctions;
  set<string> _notActiveFunctions;

  map<string, int> _clusterEdgesToNumActiveFunctions;
  map<string, int> _notClusterEdgesToNumActiveFunctions;

  int _numClusterEdgesWithSomeActiveFunction;
  int _numNotClusterEdgesWithSomeActiveFunction;

  void _initialize();

  void _activateFunction(string func);

  void _deactivateFunction(string func);

  long double _calcGroupHypPvalue() const;

public:
  EdgeGroupHyp(const MyGraph& clusterGraph, const MyGraph& universeGraph,
               MyAnnotations& annotations, const GeneOntology& go, string category);

  virtual ~EdgeGroupHyp()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  // not used
  virtual void computeEnrichmentForFunction(const string funcType, const string funcId)
    {}

  virtual void printEnrichmentResults(ostream &ostr);

};





class FuncAssociate: public FunctionalEnrichmentAlgorithm
{
private:

  MyPointSet _clusterRanks;

  vector< pair< MyNodeId, MyNT > > _currClusterRank;
  string _currClusterId;

  ofstream _outputStream;

  struct FunctionInfo{
    BioFunction function;
    set< MyNodeId > funcGenesInUniverse;
    unsigned int numFuncGenesInCluster;
    double hypPval;
    double bestHypPval;
    set< MyNodeId > genesInBestPrefix;
    set< MyNodeId > funcGenesInBestPrefix;
    pair< unsigned int, unsigned int > adjustedPval;
  };

  void computeBestPval(const BioFunction* function);

  void computeAdjustedPvals();

  // Updates _currClusterRank with the ranking given by sample i in _clusterRanks
  void updateCurrClusterInfo(unsigned int i);

  void printHeader(string sample);

  void printEnrichmentForFuncInfo(FunctionInfo &fi);

  /*
  void printEnrichmentForSample(ostream &ostr, string sample);
  */

public:
  FuncAssociate(const MyPointSet &clusterPoints, const set<string> &universeSet,
                MyAnnotations& annotations, const GeneOntology& go, string category,
                string outfile);

  virtual ~FuncAssociate()
    {}

  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  // not used
  virtual void computeEnrichmentForFunction(const string funcType, const string funcId)
   {}

  virtual void printEnrichmentResults(ostream &ostr);

};




class SetCoverEnrichment: public FunctionalEnrichmentAlgorithm
{
private:

  MyGraph _universeGraph;
  MyGraph _clusterGraph;

  set<string> _currClusterEdges;

  struct FuncStats{
    set<string> inEdges;
    set<string> outEdges;
    MyNT score;
  };

  map<string, FuncStats> _stats;

  vector<string> _selectedFunctions;
  set<string> _selectedInEdges;
  set<string> _selectedOutEdges;



  void _initialize();

  MyNT _score(unsigned int in, unsigned int out);

  string _getHighestScoringFunc();

  void _removeFunc(string func);


public:
  SetCoverEnrichment(const MyGraph& clusterGraph, const MyGraph& universeGraph,
                 MyAnnotations& annotations, const GeneOntology& go, string category);

  virtual ~SetCoverEnrichment()
    {}


  virtual void computeEnrichedFunctions(const set< string >& functionsToProcess);

  virtual void computeEnrichmentForFunction(const string funcType, const string funcId)
    { }

  virtual void printEnrichmentResults(ostream &ostr);

};







class FunctionalEnrichmentAlgorithmDispatcher
{
public:

  enum FunctionalEnrichmentAlgorithmType
    {
      NO_ALGORITHM,
      FUNCASSOCIATE,
      GENGO,
      GROUP_HYP,
      EDGE_GROUP_HYP,
      HYPERGEOMETRIC,
      PAGE,
      RANDOM_FUNCTION,
      RANDOM_UNIVERSE,
      SET_COVER,
      MAX_ALGORITHM
    };

  static FunctionalEnrichmentAlgorithm* getFunctionalEnrichmentAlgorithm(
        unsigned int algorithmType,
        const set<MyNodeId> &clusterSet, const MyGraph &clusterNet, const MyPointSet &clusterPoints,
        const set<MyNodeId> &universeSet, const MyGraph &universeNet, const MyPointSet &universePoints,
        MyAnnotations &annotations, const GeneOntology &go, string category,
        unsigned int numRandomNets, unsigned int numRandomizationIterations,
        double genGOp, double genGOq, double genGOpenalty, string outfile);

};


class EnrichmentPostprocess
{
private:

  struct enrichmentResult{
    GOFunction* function;
    double lexComponentsByNodesPValue;
    double lexComponentsByEdgesPValue;
    double densityPValue;
  };

  GeneOntology * _go;

  vector< GOFunction* > _topoSortedGOFunctions;

  double _pvalCutoff;

  map< string, enrichmentResult > _results;

  map< string, set< string > > _functionsBelowCutoff;


  map< string, set< string > > _leavesBelowCutoff;

public:
  EnrichmentPostprocess(GeneOntology& go, string resultsFile, double pvalCutoff);

  ~EnrichmentPostprocess()
    {}

  void readEnrichment(string inFile);

  void print(ostream &ostr);

};


#endif // _ENRICHMENT_ALGORITHM_H