netman.h

/**************************************************************************
 * Copyright (c) 2005-2011 T. M. Murali                                   *
 * Copyright (c) 2009-2011 Christopher L. Poirel                          *
 * Copyright (c) 2011 Christopher D. Lasher                               *
 * Copyright (c) 2010 Jacqueline Addesa                                   *
 * Copyright (c) 2003-2005 Deept Kumar                                    *
 *                                                                        *
 * 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 _NETMAN_H_
#define _NETMAN_H_


// for transform in NetMan::getCorrelationMethod().
#include <algorithm>
#include <set>
#include <string>

#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <math.h>

using namespace std;

#include "boost/multi_index_container.hpp"
#include "boost/multi_index/member.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 MyAnnotations.
#include "old-annotations.h"
#include "enrichment.h"
// for GeneOntology
#include "GO.h"
#include "histogram.h"

#include "active-networks.h"
#include "active-networks-opts.h"

// add '/' at the end, if necessary.
inline string _normaliseDirectory(string dir)
{
  if ('/' != dir[dir.length() - 1])
    dir += '/';
  return(dir);
}


void _reverseAliases(const map< string, set< string > > &nodeAliases,
                     map< string, set< string> > &reversedAliases);



struct PPNodeStats
{
public:
  string _nodeId;
  MyNT _tstat;
  MyNT _pvalue;
  MyNT _zscore;
public:
  PPNodeStats (string n, MyNT t, MyNT p, MyNT z)
      : _nodeId(n), _tstat(t), _pvalue(p), _zscore(z)
    {}
};

struct pa_nodeid_index {};
struct pa_zscore_index {};

typedef multi_index_container<
  PPNodeStats,
  indexed_by<

    // define single-column primary key.
    ordered_unique<
    tag< pa_nodeid_index >,
    member< PPNodeStats, string, &PPNodeStats::_nodeId > >,

    // there can be many nodes with the same z-score.
    ordered_non_unique<
    tag< pa_zscore_index >,
    member< PPNodeStats, MyNT, &PPNodeStats::_zscore > >

    > // matches indexed_by
  > SetOfPPNodeStats;

class NetMan
{
private:

  gengetopt_args_info _options;
  string _commandLine;

  // -a
//  vector< string > _activeNetworkFiles;

  // -b, --compute-betweenness
  bool _computeBetweennessCentralities;

  // -C
  MyNT _correlationThreshold;
  // currently, i set this value based on the value of
  // _correlationThreshold in setCorrelationThreshold().
  bool _useRandomisedGEData;


  // -E
  string _activeNetworkEdgesFile;

  // --edge-annotations-file
  MyAnnotations _edgeAnnotations;
  string _edgeAnnotationsFile;


  // -f, --functions-file
  MyAnnotations _nodeAnnotations;
//  string _nodeAnnotationsFile;
  vector< string > _nodeAnnotationsFiles;


  MyGraph wiringDiagram;
  bool _computeDenseFunctions;
  map< string, MyGraph > _denseFunctionSubgraphs;
  MyGraph _allDenseFunctionSubgraphs;
  bool _useMaximumSpanningTree;
  bool _deleteDenseFunctionSubgraphsCompletely;

  map< string, MyGraph > functionSubgraphs;

  // store all ge data in this variable. 2008-12-12, start the
  // transition from vector< GEData >.
  GEDataSets geDataSets;
public:
  vector< GEData > geData;
private:
  // map storing which class is in which index.
  map< string, unsigned int > _classToIndex;

  string currentExperimentName;
  vector< GEData >::iterator currentGEData;
  vector< string > directories;
  vector< string > experimentNames;

  // -I
  set< string > _ignoredEdgeTypes;
  string _ignoredEdgeString;
  MyGraph _ignoredWiringDiagram;

  // -l
  bool _computeNetworkLegos;


  // -m, --multiple-hypotheses-correction
  LIBENRICHMENT_TEST_TYPE _multipleHypothesesCorrectionTest;

  // --min-number-conditions
  unsigned int _minNumberConditions;

  // --min-number-interactions
  unsigned int _minNumberInteractions;

  // -N
  string _activeNetworkNodesFile;

public:
  // --node-alias-file
  map< string, set< string > > _nodeAliases;
  map< string, set< string > > _reversedNodeAliases;
private:
  
  // --node-url
  string _nodeURL;

  // -r, --num-random-runs
  unsigned int _numRandomRuns;
  bool _runOnlyRandomRuns;

  // -s
  bool _splitActiveNetworksIntoComponents;

  // --translate-annotations.
  bool _translateAnnotations;

  // --www
  string _currentWWWDirectory;

  // file containing information on the class each chip in a gene
  // expression data set belongs to.
  string classFile;
  string functionSubgraphsFile;
  string geFile;
  string networkFile;
  string oFile;

  GeneOntology _goDAG;
  bool _goDAGRead;
  string _oboFile;

  // -p
  MyNT _pvalueThreshold;

  // -t
  MyNT _edgeWeightThreshold;

  // -T
  bool _useStudentsTDistribution;

  // -x
  bool _expandActiveNetworks;


  // -z
  bool _computeActiveNetworksTillTheBitterEnd;

public:
  
  // each of these output files will contain tab-delimited data about
  // the ActiveNetworks in a form fit for loading into a database.


  // to output in 0-1 matrix format.
  ofstream binaryMatrixStream;
  // to output condition in each ActiveNetwork.
  ofstream conditionStream;

  // the edges in an ActiveNetwork or in a combined ActiveNetwork.
  ofstream edgeStream;
  // statistics on each edge in an ActiveNetwork or in a combined ActiveNetwork.
  ofstream edgeCentralityStream;

  // statistics on functional enrichment for edges.
  ofstream edgeEnrichmentStream;
  // for each edges, which functions annotating it are enriched.
  ofstream edgeEnrichedFunctionStream;
  // graph containing edge similarities of ANs.
  ofstream edgeSimilarityGraphStream;
  // the count of edge types in an ActiveNetwork or in a combined ActiveNetwork.
  ofstream edgeTypesCountStream;
  // histogram of edge weights in ANs/NLs.
  ofstream edgeWeightsHistogramStream;


  // statistics on function subgraphs.
  ofstream functionSubgraphsStream;

  // histogram info.

  // the histogram of expression values in a gene expression dataset.
  ofstream expressionValuesHistogramStream;
  ofstream histogramStream;

  // lattice edges.
  ofstream latticeStream;
  // transitively closed lattice.
  ofstream closedLatticeStream;

  // generic information.
  ofstream logStream;
  // the nodes in an ActiveNetwork or in a combined ActiveNetwork.
  ofstream nodeStream;
  // statistics on each node in an ActiveNetwork or in a combined ActiveNetwork.
  ofstream nodeCentralityStream;

  // statistics on functional enrichment for nodes.
  ofstream nodeEnrichmentStream;
  ofstream nodeEnrichmentOntologizerStream;
  // for each node, which functions annotating it are enriched.
  ofstream nodeEnrichedFunctionStream;
  // analysis of functional enrichment of nodes across ANs.
  ofstream processedNodeEnrichmentsStream;
  ofstream processedNodeEnrichmentsOntologizerStream;
  // graph containing node similarities of ANs.
  ofstream nodeSimilarityGraphStream;
  // histogram of node weights in ANs/NLs.
  ofstream nodeWeightsHistogramStream;

  // to output in tab-delimited "itm" format.
  ofstream itemsetStream;
  // statistics about ActiveNetworks in random graphs.
  ofstream randomStream;
  ofstream randomPvaluesStream;
  // various stats about ActiveNetworks.
  ofstream statsStream;

  // string streams to print info to before opening output files.
  stringstream logStringStream;
//  stringstream expressionValuesHistogramStringStream;
//  stringstream expressionValuesGnuplotStringStream;

  // the names of the output files. the first string is the type of
  // the file (same as the string in outputStreams. the second string
  // is the actual file name.
  map< string, string > _outputFileNames;
  // a convenience variable to simplify the code for opening and
  // closing all the above streams.
  map< string, ofstream* > outputStreams;

private:
  set< string > _negativeWeightEdgeTypes;

  MyHistogram _expressionValuesHistogram;

  // am i re-processing previously computed ANs?
  bool _reprocessComputedActiveNetworks;
  // am i re-processing previously computed NLs?
  bool _reprocessComputedNetworkLegos;

  // map from treatment to set of controls. each treatment can have
  // >=1 control, but each t-c pair is analysed separately. this data
  // structure is also used in Netman::createWebPages() (and therefore
  // in SetOfActiveNetworks::layout() to decide which conditions to
  // group together when drawing function-specific networks.
  map< string, set< string > > _treatmentControlPairs;


public:

  NetMan(void);
  ~NetMan(void);

  /************ START OF COMMAND LINE RELATED METHODS. *************/

  // -a
  // void addActiveNetworkFile(string optarg)
  //   {
  //     _activeNetworkFiles.push_back(optarg);
  //   }
  // unsigned int getNumActiveNetworkFiles() cons
  //   {
  //     return(_activeNetworkFiles.size());
  //   }

  // string getActiveNetworkFiles(unsigned int index) const
  //   {
  //     return(_activeNetworkFiles[index]);
  //   }


  // void getActiveNetworkFiles(vector< string > &files) const
  //   {
  //     files = _activeNetworkFiles;
  //   }

  // MORON! calling this method set() confuses the poor compilers when
  // they see set used as an STL data structure later.

  virtual void setParameters(int argc, char **argv);

  // -b, --compute-betweenness
  void setComputeBetweennessCentralities()
    {
      _computeBetweennessCentralities = true;
    }
  bool getComputeBetweennessCentralities() const
    {
      return(_computeBetweennessCentralities);
    }



  // -B
  void setOBOFile(char* file)
    {
      _oboFile = file;
    }

  // -c
  // call this method only after setting the experiment name (NetMan::expt_name).
  void setClassFile(string name)
    {
      classFile = name;
    }
  string getClassFile(void) const
    {
      return classFile;
    }


  // -C
  void setCorrelationThreshold(float arg)
    {
      _correlationThreshold = arg;//atof(arg.c_str());
      if (1 < _correlationThreshold)
        setUseRandomisedGEData();
    }

  // --compute-stability.
  bool getComputeNetworkLegoStability() const
    {
      return(_options.compute_stability_given);
    }

  void setUseRandomisedGEData()
    {
      _useRandomisedGEData = true;
    }
  bool getUseRandomisedGEData() const
    {
      return(_useRandomisedGEData);
    }

  // --correlation-method
  string getCorrelationMethod() const
    {
      string method = _options.correlation_method_arg;
      // see http://www.devx.com/getHelpOn/Article/9702/1954?pf=true
      transform(method.begin(), method.end(), method.begin(), (int(*)(int)) tolower);
      return(method);
//      return(_options.correlation_method_arg);
    }

  MyNT getCorrelationThreshold() const
    {
      return(_correlationThreshold);
    }


  // -d
  void addDirectory(char* dir)
    {
      directories.push_back(dir);
    }
  unsigned int getNumDirectories() const
    {
      return(directories.size());
    }
  string getDirectory(unsigned int index) const
    {
      return(directories[index]);
    }



  // -D
  void setComputeDenseFunctions()
    {
      _computeDenseFunctions = true;
    }
  bool getComputeDenseFunctions() const
    {
      return(_computeDenseFunctions);
    }
  
  // -E
  void setActiveNetworksEdgesFile(string file)
    {
      _activeNetworkEdgesFile = file;
    }
  string getActiveNetworksEdgesFile() const
    {
      return(_activeNetworkEdgesFile);
    }

  bool getProcessComputedActiveNetworks() const
    {
      return(("" != getActiveNetworksEdgesFile())
             || ("" != getActiveNetworksNodesFile()));;
    }

  // --experiment-name
  string getExperimentName() const
    {
      return(_options.experiment_name_arg);
    }

  // -f
  void setNodeAnnotationsFile(char* name)
    {
//      _nodeAnnotationsFile = name;
      _nodeAnnotationsFiles.push_back(name);
    }

  // -g
  void setGEfile(char* name)
    {
      geFile = name;
    }
  string getGEfile(void) const
    {
      return geFile;
    }


  // -I
  void addIgnoredEdgeType(string type)
    {
      // _ignoredEdgeString must start with a "-" if it is non-empty.
      if (0 == _ignoredEdgeTypes.size())
        _ignoredEdgeString = "-NO";
      _ignoredEdgeString += "-" + type;
      _ignoredEdgeTypes.insert(type);
    }

  //  am I computing ANs?
  bool getComputeActiveNetworks(void) const
    {
      //  I am computing ANs if ..
      //
      // i am not computing NLs
      return(!getComputeNetworkLegos()
              // and i am not computing pathway activation
              && !getComputePathwayPerturbation()
              // and i am not just processing the network.
              && !getProcessNetwork()
              // and i don't have a pre-computed AN file.
              && !getProcessComputedActiveNetworks());
    }


  // -l
  void setComputeNetworkLegos()
    {
      _computeNetworkLegos = true;
    }

  bool getComputeNetworkLegos() const
    {
      return(_computeNetworkLegos);
    }

  // -L
  bool setFunctionSubgraphsFile(string arg)
    {
      functionSubgraphsFile = arg;
    }

  // -m, --multiple-hypotheses-correction
  void setMultipleHypothesesCorrectionTest(string optarg)
    {
      _multipleHypothesesCorrectionTest =
        convertStringToTestType(optarg);
    }
  LIBENRICHMENT_TEST_TYPE getMultipleHypothesesCorrectionTest() const
    {
      return(_multipleHypothesesCorrectionTest);
    }


  
  // --min-number-conditions
  void setMinNumberConditions(unsigned int min)
    {
      _minNumberConditions = min;
    }
  unsigned int getMinNumberConditions() const
    {
      return(_minNumberConditions);
    }

  // --min-number-interactions
  void setMinNumberInteractions(unsigned int min)
    {
      _minNumberInteractions = min;
    }
  unsigned int getMinNumberInteractions() const
    {
      return(_minNumberInteractions);
    }

  // -n
  void setNetworkFile(char* file)
    {
      networkFile = file;
    }

  string getNetworkFile(void) const
    {
      return networkFile;
    }

  // -N
  void setActiveNetworksNodesFile(string file)
    {
      _activeNetworkNodesFile = file;
    }
  string getActiveNetworksNodesFile() const
    {
      return(_activeNetworkNodesFile);
    }

  bool getProcessNetwork() const
    {
      return(//(0 == getNumActiveNetworkFiles())
        //&&
             (0 == getNumDirectories())
             && ("" == getGEfile())
             && ("" != getNetworkFile())
             );
    }


  // --layout-constraints-file
  string getLayoutConstraintsFile() const
    {
      return(_options.layout_constraints_file_arg);
    }


  // -p
  void setPvalueThreshold(float thresh)
    {
      _pvalueThreshold = thresh;//atof(thresh.c_str());
    }
  MyNT getPvalueThreshold() const
    {
      return(_pvalueThreshold);
    }

  // --pathway-perturbation, -P
  bool getComputePathwayPerturbation() const
    {
      return(_options.pathway_perturbation_given);
    }
      
  // -r
  void setNumRandomRuns(unsigned int num)
    {
      _numRandomRuns = num;//atoi(arg);
    }


  void setOnlyRandomRuns()
    {
      _runOnlyRandomRuns = true;
    }

  bool getOnlyRandomRuns() const
    {
      return(_runOnlyRandomRuns);
    }

  unsigned int getNumRandomRuns() const
    {
      return(_numRandomRuns);
    }


  // --regulation-bound.
  MyNT getRegulationBound() const
    {
      if (_options.regulation_bound_given)
        return(_options.regulation_bound_arg);
      return(1);
    }
  MyNT getUpRegulationBound() const
    {
      return(getRegulationBound());
    }
  MyNT getDownRegulationBound() const
    {
      return(-getRegulationBound());
    }

  string getRegulationMethod() const
    {
      return(_options.regulation_method_arg);
    }

  // --sa-factor-number-iterations.
  unsigned int getSimulatedAnnealingFactorNumberIterations() const
    {
      return(_options.sa_factor_number_iterations_arg);
    }


  // --sparsify-network
  bool getSparsifyNetwork() const
    {
      return(_options.sparsify_network_given);
    }
  string getSparsificationMethod() const
    {
      return(_options.sparsify_network_arg);
    }

  // -s, --split-active-networks
  void setSplitActiveNetworksIntoComponents()
    {
      _splitActiveNetworksIntoComponents = true;
    }
  bool getSplitActiveNetworksIntoComponents() const
    {
      return(_splitActiveNetworksIntoComponents);
    }

  // --split-gene-expression-data
  bool getSplitGEData() const
    {
      return(_options.split_gene_expression_data_given);
    }

  // -t
  void setEdgeWeightThreshold(string optarg)
    {
      _edgeWeightThreshold = atof(optarg.c_str());
    }
  void setEdgeWeightThreshold(MyNT thresh)
    {
      _edgeWeightThreshold = thresh;
    }
  MyNT getEdgeWeightThreshold() const
    {
      return(_edgeWeightThreshold);
    }


  // -T
  void setUseStudentsTDistribution()
    {
      _useStudentsTDistribution = true;
    }
  bool getUseStudentsTDistribution() const
    {
      return(_useStudentsTDistribution);
    }

  // --translateAnnotations
  bool getTranslateAnnotations() const
    {
      return(_options.translate_annotations_given);
    }

  // -trivial-network-legos
  bool getComputeTrivialNetworkLegos() const
    {
      return(_options.trivial_network_legos_given);
    }




  // -x
  void setExpandActiveNetworks()
    {
      _expandActiveNetworks = true;
    }
  bool getExpandActiveNetworks() const
    {
      return(_expandActiveNetworks);
    }

  // -z
  void setComputeActiveNetworksTillTheBitterEnd()
    {
      _computeActiveNetworksTillTheBitterEnd = true;
    }
  bool getComputeActiveNetworksTillTheBitterEnd() const
    {
      return(_computeActiveNetworksTillTheBitterEnd);
    }


  // WWW-related options.

  // --embed
  bool getEmbedActiveNetworkInSuperGraph() const
    {
      return(_options.embed_given);
    }

  // --www
  bool getCreateWebPages() const
    {
#ifdef CTEMPLATE
      return(_options.www_given);
#else
      return(false);
#endif
    }

  void createWebPages(SetOfActiveNetworks &anetSet);


  // --www-params-dir
  string getWWWParametersDir() const
    {
      return(_normaliseDirectory(_options.www_params_dir_arg));
    }

  /************ END OF COMMAND LINE RELATED METHODS. *************/

  void getWiringDiagram(MyGraph &wd) //const: COMPILATION ERROR WITH const.
    {
      wd = wiringDiagram;
    }

  void getWiringDiagramNodes(set< string > &nodes) const
    {
      wiringDiagram.getNodeSet(nodes);
    }

  string getEdgeAnnotationsFile(void) const
    {
      return(_options.edge_annotations_file_arg);
    }

//  string getNodeAnnotationsFile(void) const;
  void getNodeAnnotationsFiles(vector< string > &files) const;
  string getNodeAnnotationsFile(unsigned int index) const
    {
      return(_nodeAnnotationsFiles[index]);
    }
  unsigned int getNumNodeAnnotationsFiles() const
    {
      return(_nodeAnnotationsFiles.size());
    }

  // --node-alias-file
  string getNodeAliasFile() const
    {
      return(_options.node_alias_file_arg);
    }


  string getFunctionSubgraphsFile() const
    {
      return(functionSubgraphsFile);
    }

  string getOBOFile(void) const
    {
      return(_oboFile);
    }

  // --treatment-control-file
  string getTreatmentControlFile() const
    {
      return(_options.treatment_control_file_arg);
    }


  void readEdgeAnnotations(string name);

  void readNodeAliases(string file);

  void readNodeAnnotations(string name, set< string > *geneSet = NULL);

  void readOBO(string name)
    {
      _goDAG.read(name);
      // i better do this if i want depths to be printed out.
      _goDAG.computeFunctionDepths();
      _goDAGRead = true;
    }

  bool readFunctionSubgraphs(string);
  bool readGEData(string);
  bool readNetwork(string, MyNT minimumEdgeWeight = 0);

  vector< SetOfPPNodeStats > treatAliases(MyGraph &pathway, const GEData &geData);

  bool allowNegativeEdgeWeight(set< string >& types) const
    {
      set< string >::const_iterator itr;
      for (itr = _negativeWeightEdgeTypes.begin(); itr != _negativeWeightEdgeTypes.end();
           itr++)
        if (types.end() != types.find(*itr))
          return(true);
      return(false);
    }

  // void combineActiveNetworks(vector< string > &activeNetworkFiles);



  void computeDenseFunctions(ostream *fstr,
                             map< string, MyGraph > &denseFunctionSubgraphs);

  void computeDenseFunctions(ostream *fstr = NULL)
    {
      computeDenseFunctions(fstr, _denseFunctionSubgraphs);
      map< string, MyGraph >::iterator gitr;
      for (gitr = _denseFunctionSubgraphs.begin(); gitr != _denseFunctionSubgraphs.end();
           gitr++)
#ifdef DEBUG
        {
          cout << "Function with dense subgraph is " << gitr->first
               << " with " << gitr->second.numNodes() << " nodes and "
               << gitr->second.numEdges() << " edges." << endl;
          gitr->second.printEdges(cout);
#endif
          _allDenseFunctionSubgraphs.add(gitr->second);
#ifdef DEBUG
        }
      cout << "Union of all functions with dense subgraphs has "
           << _allDenseFunctionSubgraphs.numNodes() << " nodes and "
           << _allDenseFunctionSubgraphs.numEdges() << " edges." << endl;
      _allDenseFunctionSubgraphs.printNodes(cout);
      _allDenseFunctionSubgraphs.printEdges(cout);
#endif
    }

  void computePathwayPerturbation();

  // i cannot make pathways const because MyGraph has no const node iterator.
  void computePathwayPerturbation(const map< string, set< string > > &treatmentControlPairs, map< string, MyGraph> &pathways);

  void computePathwayPerturbation(MyGraph &pathway, const GEData &treatment,
                                  const GEData &control, string contrast, ActiveNetwork &subgraph);

  void computePathwayPerturbation(MyGraph &pathway, const GEData &geData);


//  static bool compareFunction(MyNode x, MyNode y);


  void plotExpressionValuesHistogram();

  void plotExpressionCorrelationsHistogram();

  void processNetworkLegos(SetOfActiveNetworks &anetSet);


  void readDirectories();

  void removeUnresponsiveGenes(GEData &gedata);

  bool computeCorrelations(MyGraph &wiring);//GEData &expression);

  // \param[in] randomRun, if true, the method knows that the wiring
  // diagram is a randomly generated network. In this case, the method
  // does not print out some statistics about the dense subgraphs it
  // finds.
  bool computeDenseSubgraphs(vector< ActiveNetwork > &denseSubgraphs, bool randomRun = false);


  void setWiringDiagramEdgeWeightsFS()
    {
      wiringDiagram.computeEdgeFSWeights();
    }

  void setWiringDiagramEdgeWeightsCD()
    {
      wiringDiagram.computeEdgeCDWeights();
    }

  void postProcess(SetOfActiveNetworks &anetSet);

  void process();

  void processComputedActiveNetworks();


  void processFunctionSubgraphs();

  void processNetwork();


  void randomise(MyHistogram &histogram);

  void separateDenseFunctions();



private:

  // assign edge weights based on node stats.
  void _assignEdgeWeights(MyGraph &pathway, SetOfPPNodeStats &nodeStats,
                          ostream &ostr);

public:
  void _openOutputFiles();
  void _closeOutputFiles();

private:
  // check if ids are same for different datasets.
  void _performSanityChecks();

  void _processDenseFunctionSubgraphs(MyGraph &wiringDiagram, MyGraph &mst);

  void _readTreatmentControlPairs(string file);

};

#endif