[Bio] / FluxMapMaker / FluxMapMakerExe.cc Repository:
ViewVC logotype

View of /FluxMapMaker/FluxMapMakerExe.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (download) (as text) (annotate)
Wed Jan 28 19:06:24 2004 UTC (15 years, 8 months ago) by efrank
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +21 -13 lines
who knows at this point

//************************************************************************
//  Generate Dot/Dotty rep of a network for graphing.
//
//
//************************************************************************

#include <stdio.h>
#include <string>
#include <vector>
#include <set>
#include <map>
#include <iterator>
#include <fstream>
#include <sys/types.h>
#include <sys/stat.h>

#include "RecoReadout/RecoElem.hh"
#include "RecoReadout/Reaction.hh"

#include "RecoReadoutOracle/BssFactory.hh"
#include "RecoReadoutOracle/OracleMgr.hh"
#include "RecoElemSBW/SbmlFactory.hh"

using namespace std;

//************************************************************************
void usage( void );
void getMultiNodeMap( map<const string,bool>& m );
void fakeGatherRxns( RecoElem& topNode, std::set<Reaction>& rxns );
void gatherRxns( RecoElem& topNode, std::set<Reaction>& rxns );
void makeFluxMap( const set<Reaction>& rxns, 
		  const map<const string,bool>& multiNodeMap );
void makeRBreaker( const set<Reaction>& rxns );

//************************************************************************
int
main(int argc, char* argv[] ) {
//************************************************************************

  std::set<Reaction>          rxns;           // collect reactions here
  std::map<const string,bool>  multiNodeMap;   // names of nodes to duplicate.
                                              //    see getMultiNodeMap()
  struct stat sb;

  if ( 2 != argc && 3 != argc) {
    usage();
    return 1;
  }

  RecoElem*  recoTop(0);
  if (0==stat( argv[1], &sb ) && S_ISREG( sb.st_mode ) ) {
    recoTop = SbmlFactory( argv[1] );
  } else {
    // Get the reconstruction out of the Database
    OracleMgr::mgr().startJob();
    BssFactory factory;
    unsigned int topNodeId = atoi( argv[1] );
    recoTop = factory.getRecoElemById( topNodeId );
    OracleMgr::mgr().endJob();
  }
      
  getMultiNodeMap( multiNodeMap );
  gatherRxns( *recoTop, rxns );
      
  if(argc == 3 && strcmp(argv[2], "-dot") == 0)
	  makeFluxMap( rxns, multiNodeMap );
  else
	  makeRBreaker( rxns );

  delete recoTop;
  OracleMgr::mgr().endJob();

  return 0;
}

//************************************************************************
void
usage( void ) {
//************************************************************************
  cout << endl;
  cout << "usage: FluxMapMaker TopNodeId [-dot]" << endl;
  cout << "   TopNodeId is the ProcessId of the top of the reconstruction." << endl;
  cout << "   You can find the top node of a reconstruction with lsCat," << endl;
  cout << "   which lists the Catalog." << endl;
  cout << "   If -dot is specified, FluxMapMaker will output directly to " << endl;
  cout <<"    a .DOT file instead of an Rbreaker file." << endl;
  cout << endl;

  return;
}



//************************************************************************
void
fakeGatherRxns( RecoElem& topNode,
	    std::set<Reaction>& rxns ) {
//************************************************************************
  
  // fake a test set
  //    A + B -> C + D
  //    D + E -> P

  Reaction  r1( 1, Reaction::forward);
  Reaction r2( 2, Reaction::both);

  vector<string>& r1In =  r1.getInputMetaboliteNames();
  vector<string>& r1Out=  r1.getOutputMetaboliteNames();
  vector<int>& r1InStoich =  r1.getInputStoich();
  vector<int>& r1OutStoich =  r1.getOutputStoich();

  vector<string>& r2In =  r2.getInputMetaboliteNames();
  vector<string>& r2Out = r2.getOutputMetaboliteNames();
  vector<int>& r2InStoich =  r2.getInputStoich();
  vector<int>& r2OutStoich =  r2.getOutputStoich();

  r1In.push_back( "A" );
  r1InStoich.push_back( 1 );
  r1In.push_back( "B" );
  r1InStoich.push_back( 1 );
  r1Out.push_back( "C" );
  r1OutStoich.push_back( 1 );
  r1Out.push_back( "D" );
  r1OutStoich.push_back( 1 );


  r2In.push_back( "D" );
  r2InStoich.push_back( 1 );
  r2In.push_back( "E" );
  r2InStoich.push_back( 1 );

  r2Out.push_back( "P" );
  r2OutStoich.push_back( 1 );

  rxns.insert( r1 );
  rxns.insert( r2 );

  cout <<  "fake1: " << r1.asString( ) << endl;
  cout <<  "fake2: " << r2.asString( ) << endl;

  return;
}

//************************************************************************
void
gatherRxns( RecoElem& topNode,
	    std::set<Reaction>& rxns ) {
//************************************************************************
  

  std::vector<RecoElem>& kids    = topNode.getSubElements();
  int                    i;
  int                    nKids   = kids.size();
  std::vector<Reaction>& myRxns  = topNode.getReactions();
  int                    nRxn    = myRxns.size();
  
  // recurse to gather rxns from children
  for( i=0; i<nKids; i++ ) gatherRxns( kids[i],  rxns );
  

  // unwind by puttin my rxns onto the set.
  
  for( i=0; i<nRxn; i++) rxns.insert( myRxns[i] );

  return;
}

//************************************************************************
void
getMultiNodeMap( map<const string,bool>& m ) {
//************************************************************************
  // Some nodes in the reaction network are special, e.g., ATP is so
  // ubiquitous that you do not want just one ATP node that ties to all
  // the places it is used.  The result would be a tangle.  Instead, some
  // nodes are replicated every time they are used...like a call-out or
  // reference.
  //
  // Get the list of these.

  // really, this should come from a file but for now we hack it in to
  // get started.


  //ADP(red)
  //ATP(green)

  m.insert( pair<string,bool>( "ADP", true ));
  m.insert( pair<string,bool>( "ATP", true ));
  m.insert( pair<string,bool>( "AMP", true ));
  m.insert( pair<string,bool>( "GTP", true ));
  m.insert( pair<string,bool>( "GDP", true ));
  m.insert( pair<string,bool>( "GMP", true ));
  m.insert( pair<string,bool>( "H2O", true ));
  m.insert( pair<string,bool>( "H2O2", true ));
  m.insert( pair<string,bool>( "O2", true ));
  m.insert( pair<string,bool>( "PI", true ));
  m.insert( pair<string,bool>( "COA", true ));
  //m.insert( pair<string,bool>( "ACCOA", true ));
  m.insert( pair<string,bool>( "CO2", true ));
  m.insert( pair<string,bool>( "NAD", true ));
  m.insert( pair<string,bool>( "NADP", true ));
  m.insert( pair<string,bool>( "NADH", true ));
  m.insert( pair<string,bool>( "NADPH", true ));
  m.insert( pair<string,bool>( "QH2", true ));
  m.insert( pair<string,bool>( "Q", true ));
  //m.insert( pair<string,bool>( "PYR", true ));
  m.insert( pair<string,bool>( "NH3", true ));


  return;
}


//************************************************************************
void
makeRBreaker( const set<Reaction>& rxns ) {
//************************************************************************

  std::set<Reaction>::iterator iRx		= rxns.begin();
  std::set<Reaction>::iterator iEnd		= rxns.end();
  vector<string>               enzymeNames;
  std::string                  printEnzymeName;

  // this stuff is the RBreaker magic mantra for "duplicate these
  // nodes"  Notice that this assumes specific names.  Ramya is
  // fixing with her DB.

  cout << "X\tADP(red):ATP(green):AMP:GTP:GDP:GMP:H2O:H2O2:O2" << endl
       << "X\tPI:COA:ACCOA:CO2" << endl
       << "X\tNAD:NADP:NADH:NADPH:QH2:Q:PYR:NH3" << endl;

  for( ; iRx != iEnd; iRx ++ ) {
    enzymeNames = iRx->getEnzymeNames();
    cout << "R\t"
	 << "R" << iRx->getId() << "\t"
	 << "geneId\t"
	 << iRx->asString() << "\t";

//    cout << "Rx" << iRx->getId() << "\t";
    cout << iRx->asString() << "\t";
	 
    // if ( 0 == enzymeNames.size() ) {
    //   cout << "Rx" << iRx->getId() << "\t";
    // } else {
    //  cout << enzymeNames[0] << ":Rx" << iRx->getId() << "\t";
    // }
    cout << endl;
  }

  return;
}



//************************************************************************
void
makeFluxMap(const set<Reaction>& rx, const map<const string,bool>& multiNodeMap ) {
//************************************************************************

  ofstream theFile;
  theFile.open( "fluxmap.dot", ofstream::trunc);

  // header in dot file

  theFile << "digraph PW {" << endl;
  theFile << "graph [rankdir=LR];" << endl;
  theFile << "node [shape=box, height=0.1]" << endl;


  // the reactions to process

  std::set<Reaction>::iterator i     = rx.begin();
  std::set<Reaction>::iterator iEnd  = rx.end();

  // the set of species for which we've made at least one node.  the int
  // is the node number we assigned.

  std::map<string,int> speciesMap;

  int  nodeId=0;
  int  thisNodeId;
  bool multiNode;
  bool createNewNode;
  
  for( ; i != iEnd; i++ ) {
    const vector<string>& inney = i->getInputMetaboliteNames();
    const vector<string>& outey = i->getOutputMetaboliteNames();

    // for each rx, we make a bogus pair of nodes.  All inputs point to the first
    // invisible node.  All outputs flow from the 2nd invis node.  The two are
    // connected by an edge we can label.

    int inNode = nodeId++;
    int outNode= nodeId++;

    theFile << "node" << inNode  << " [label=\"\", style=invis];" << endl;
    theFile << "node" << outNode << " [label=\"\", style=invis];" << endl;
	
    //*** input metabs

    for ( int in=0; in<inney.size(); in++ ) {
      const string& inName = inney[in];

      multiNode     = false;
      createNewNode = false;

      if ( speciesMap.find( inName ) != speciesMap.end() ) {
         if ( multiNodeMap.find( inName ) != multiNodeMap.end() ) {
	         multiNode    = true;
	          createNewNode = true;
	      }
      } else {
	      createNewNode = true;
      }

      if ( createNewNode ) {
         thisNodeId = nodeId++;
         speciesMap[ inName ] = thisNodeId;
	
         if (multiNode) {
            theFile <<  "node" << thisNodeId << " [label = \"" << inName 
            << "\"];" << endl;
         } else {
            theFile << "node" << thisNodeId << " [label = \"" << inName 
            << "\", color = \"light blue\", style=filled];" << endl;
         }
      } else {
	      thisNodeId = speciesMap[ inName ];
      }
      theFile << "node" << thisNodeId << " -> node" << inNode << ";" << endl;
    }

    //*** connecting edge

    theFile <<  "node" << inNode << " -> node" << outNode 
	    << " [label=\"" << i->asString() << "\"];" << endl;

    //*** output metabs

    for ( int out=0; out<outey.size(); out++ ) {
      const string& outName = outey[out];

      multiNode     = false;
      createNewNode = false;

      if ( speciesMap.find( outName ) != speciesMap.end() ) {
	      if ( multiNodeMap.find( outName ) != multiNodeMap.end() ) {
	         multiNode    = true;
	         createNewNode = true;
	      }
      } else {
         createNewNode = true;
      }

      if ( createNewNode ) {
	thisNodeId = nodeId++;
	speciesMap[ outName ] = thisNodeId;
	
	if (multiNode) {
	  theFile <<  "node" << thisNodeId << " [label = \"" << outName 
		  << "\"];" << endl;
	} else {
	  theFile << "node" << thisNodeId << " [label = \"" << outName 
		  << "\", color = \"light blue\", style=filled];" << endl;
	}
      } else {
	thisNodeId = speciesMap[ outName ];
      }
      theFile << "node" << outNode << " -> node" << thisNodeId << endl;
    }
  }


  // trailer in dot file

  theFile << "}" << endl;
  theFile.close();

  return;
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3