[Bio] / FluxAnalyzerGen / FluxAnalyzerGen.cc Repository:
ViewVC logotype

View of /FluxAnalyzerGen/FluxAnalyzerGen.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (download) (as text) (annotate)
Wed May 21 20:39:08 2003 UTC (17 years ago) by efrank
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +18 -11 lines
Now reads geometry file rather than assuming it is making text based maps.

//************************************************************************
// Play with exploring the e. coli reco in biosimscratch.
//
//
// look at enscript | ps2pnm | ppmtopcx to get faux fluxmap.
//
//************************************************************************

#include <stdio.h>
#include <string>
#include <vector>
#include <map>
#include <set>
#include <iterator>
#include <fstream>

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

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

using namespace std;

//************************************************************************
void usage( void );
void gatherRxns( RecoElem& topNode, std::set<Reaction>& rxns);
void gatherSpecies( const std::set<Reaction>& rx, std::set<string>& species );
void dumpFluxAnalyzerFiles( const std::set<Reaction>& rx, 
			    std::set<string>& species,
			    const std::string& mapImageFile,
			    std::map<int, pair<int, int> >& );
void  genAppPara( const std::string& mapImageFile );
void  genAssembly(const std::set<Reaction>& rx, std::set<string>& species );
void  genMacroSynthesis(const std::set<Reaction>& rx, std::set<string>& species );
void  genMacroMolecules(const std::set<Reaction>& rx, std::set<string>& species );
void  genMetabolites(const std::set<Reaction>& rx, std::set<string>& species );
void  genReactions(const std::set<Reaction>& rx, 
		   const std::set<string>& species,
		   std::map<int, pair<int,int> >& geomMap);
void  readGeomFile( const std::string& theFileName, map<int,pair<int,int> >& );


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

  std::set<string>              species;      // collect the species list here
  std::set<Reaction>            rxns;         // collect reactions here
  std::map<int, pair<int,int> > geomMap;       // rxId -> (x,y)

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

  unsigned int topNodeId = atoi( argv[1] );
  std::string  mapFileBase = argv[2];
  std::string  mapFileGeom = mapFileBase + ".geom";
  std::string  mapFileGif  = mapFileBase + ".gif";

  cout << "FluxAnalyzerGen: going for node " << topNodeId << endl;
  cout << "FluxAnalyzerGen: will use mapFile " << mapFileGeom << endl;


  // Get the reconstruction out of the Database
  OracleMgr::mgr().startJob();
  BssFactory theFactory;
  RecoElem* recoTop = theFactory.getRecoElemById( topNodeId);


  // Read the geometry
  readGeomFile( mapFileGeom, geomMap );


  // Gather the list of reactions and species

  gatherRxns( *recoTop, rxns );
  cout << "Found " << rxns.size() << " reactions." << endl;

  //  set<Reaction>::iterator i = rxns.begin();
  //  set<Reaction>::iterator e = rxns.end();
  //  for( ; i!= e; i++ ) cout << i->asString() << endl;
  

  gatherSpecies( rxns, species );
  cout << "Found " << species.size() << " metabolites." << endl;


  // Dump the FluxAnalyzer files
  
  dumpFluxAnalyzerFiles( rxns, species, mapFileGif, geomMap );


  OracleMgr::mgr().endJob();
  delete recoTop;
  return 0;
}

//************************************************************************
void
usage( void ) {
//************************************************************************
  cout << endl;
  cout << "usage: FluxAnalyzerGen TopNodeId MapFileName" << endl;
  cout << 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 << "   MapFileName is the name of the file with the metabolic map."<<endl;
  cout << "   FluxAnalyzerGen will look for MapFileName.geom and MapFileName.gif" << endl;
  cout << 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
gatherSpecies( const std::set<Reaction>& rx, std::set<string>& species ) {
//************************************************************************

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

  for(    ; i != iEnd; i++ ) {
    const std::vector<string>& inames = i->getInputMetaboliteNames();
    n = inames.size();
    for( j=0; j<n; j++) species.insert( inames[j] );
    
    const std::vector<string>& onames = i->getOutputMetaboliteNames();
    n = onames.size();
    for( j=0; j<n; j++) species.insert( onames[j] );

  }
}

//************************************************************************
void
dumpFluxAnalyzerFiles( const std::set<Reaction>& rxns,
		       std::set<string>& species,
		       const std::string& mapFileImage,
		       std::map<int, pair<int,int> >& geomMap ) {
//************************************************************************
  
  #define RXNSperPAGE 26    // hardwired #...fix later.

  unsigned int nrxns = rxns.size();
  unsigned int nMaps = nrxns / RXNSperPAGE;
  if ( 0 != nrxns % RXNSperPAGE ) nMaps++;

  genAppPara( mapFileImage );
  genAssembly( rxns, species );
  genMacroSynthesis(rxns, species );
  genMacroMolecules( rxns, species );
  genMetabolites( rxns, species );
  genReactions( rxns, species, geomMap );

}

//************************************************************************
void
genAppPara( const std::string& mapImageFile ) {
  // app_para.m has overall constants like colors and numeric limits.
  // Probably the only thing we care about is the list of flux map
  // images.
//************************************************************************
  ofstream theFile;
  theFile.open( "app_para.m", ofstream::trunc);

  std::string  widthStr;
  std::string  heightStr;
  std::string  fontStr;
  

  //  for (int i=1; i<nMaps+1; i++ ) {
  for (int i=0; i<1; i++ ) {
    widthStr += "0.08 ";
    heightStr+= "0.03 ";
    fontStr  += "8 ";
  }

  theFile << "epsilon=1e-10;" << endl
	  << "basic_color=[0.7         0.7         0.7];" << endl
	  << "cr_color=[0.5         0.5           1];" << endl
	  << "br_color=[1         0.2         0.2];" << endl
	  << "nbr_color=[0.2           1         0.2];" << endl
	  << "text_color=[0  0  0];" << endl
	  << "macro_synth_color=[1  0  0];" << endl
	  << "macro_color=[1         0.6         0.6];" << endl
	  << "box_reaction_width=[" + widthStr + "];" << endl
	  << "box_reaction_height=["+ heightStr + "];" << endl
	  << "box_macro_width=[" + widthStr + "];" << endl
	  << "box_macro_height=[" + heightStr + "];" << endl
	  << "fontsize_reaction=[" + fontStr + "];" << endl
	  << "fontsize_macro=[" + fontStr + "];" << endl
	  << "fluxmaps={" << "'FluxMap', '" << mapImageFile << "' };" << endl;

  theFile.close();
  return;
}


//************************************************************************
void
genAssembly(const std::set<Reaction>& rx, std::set<string>& species ) {
//************************************************************************
  // The assembly file associates metabolites wiht macro molecules and
  // says where to put text-boxes to show fluxes of the metabollites into
  // formation of the macromolecules.  The file does not give the
  // the stoichiometry for the processes.  That is in macromolecule_sysnthesis.
  //
  // We can not generate this file here.  It goes beyond reconstruction to
  // modeling.  The data we need are not in the reco itself.  Aha!
  ofstream theFile;
  theFile.open( "assembly", ofstream::trunc);
  theFile.close();

  return;
}

//************************************************************************
void
genMacroSynthesis(const std::set<Reaction>& rx, std::set<string>& species ) {
//************************************************************************
  // The macromolecule_synthesis file lists the stoichiometry for combining
  // metabolites into macromolecules.  The coeficients have units of
  // mmol/g and thus say how many milli-moles of each metabolite are needed
  // to make a gram (dry weight) of macromolecule.
  //
  // We can not generate this file here.  It goes beyond reconstruction to
  // modeling.  The data we need are not in the reco itself.  Aha!

  ofstream theFile;
  theFile.open( "macromolecule_synthesis", ofstream::trunc);
  theFile.close();

  return;
}

//************************************************************************
void
genMacroMolecules(const std::set<Reaction>& rx, std::set<string>& species ) {

//************************************************************************
  // The macromolecules file lists all the macromolecules giving short
  // names like "Lip" and long names like "Lipid."  The fractional composition
  // of the biomass in terms of the macromolecules is also given in this
  // file and locations of text boxes to show those values.  Note that
  // the idea here is that cell growth robs metabolites from the network
  // of reactions.  The metabolites are combined into macromolecules
  // like protein, lipid, RNA, DNA, etc.  FluxAnalyzer allows you to
  // put in phenomenological description of this by saying the relative
  // composition of the biomass in terms of these macromolecules and the
  // relative composition of the macromolecules in terms of the metabolites,
  // but you needn't list all of the biochemistry to make the macromol's.
  //
  // We can not generate this file here.  It goes beyond reconstruction to
  // modeling.  The data we need are not in the reco itself.  Aha!

  ofstream theFile;
  theFile.open( "macromolecules", ofstream::trunc);
  theFile.close();

  return;
}

//************************************************************************
void
genMetabolites(const std::set<Reaction>& rx, std::set<string>& species ) {
//************************************************************************
  // The metabolites file lists all metabolites (short name and long) in
  // the model, gives variance (???) and a flag saying whether the metabolite
  // is external (flag=1) or not.  External metabolites do not have to
  // be balanced.  

  cout << endl
       << "WARNING: The metabolites file has hardwired values for variances."
       << endl
       << "WARNING: The metabolites file has all external flags set 0 (false)"
       << endl;

  ofstream theFile;
  theFile.open( "metabolites", ofstream::trunc);

  std::set<string>::iterator i     = species.begin();
  std::set<string>::iterator iEnd  = species.end();

  for( ; i != iEnd; i++ ) {
    theFile << *i               // short name for metabolite
	    << " "		//
	    << *i 		// long name.  we don't have.  reuse short.
	    << " "		// 
	    << "0.001"		// variance.  just wing it...
	    << " "		// 
	    << "0";		// external flag. just wing it as internal.
    theFile <<endl;
  }

  theFile.close();
  return;
}

//************************************************************************
void
genReactions(const std::set<Reaction>& rxns, 
	     const std::set<string>& species,
	     std::map<int, pair<int,int> >& geomMap ) {
//************************************************************************
  // The reactions file defines the reaction network. Each line is a
  // flux.  We are supposed to give a reaction name, but generating one
  // automatically isn't trivial just given the reaction, so we just
  // use a reaction number.

  int      nMap=0;
  int      xpos=10;		// just walk boxes down page at fixed x
  int      ypos;		// compute below

  ofstream theFile;
  theFile.open( "reactions", ofstream::trunc);

  cout << endl
       << "WARNING:  Reactions file is not connected to a flux map file"
       << " and has fake information for many things." << endl;

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

  for( int nRxDone=0; i != iEnd; i++, nRxDone++ ) {

    theFile << "R" << i->getId()   // reaction name
	    << " "		   //
	    << i->asString("=").c_str() // reaction
	    << " |"		   // separator
	    << " #";	  	   // default rate.  # means "not spec'd"

    // Max an Min values for flux.
    Reaction::Direction d = i->getDirection();
    if ( d == Reaction::forward ) {
      theFile << " 0 Inf";
    } else if ( d == Reaction::backward) {
      theFile << " -Inf 0";
    } else if ( d == Reaction::both) {
      theFile << " -Inf Inf";
    }

    // get the geometry

    pair<int,int>& coord = geomMap[ i->getId() ];

    theFile << " 0"                 // coefficient for objective function
            << " " << coord.first
            << " " << coord.second
	    << " " << 1  	    // flux map number
	    << " 1"                 // editable
            << " 0.01";		    // variance???

    theFile << endl;
  }

  theFile.close();
  return;
}

//************************************************************************
void
readGeomFile( const std::string& theFileName,
	      std::map<int, pair<int,int> >& geomMap) {
//************************************************************************
  //
  // NB: geomMap< rxId, (xcoord, ycoord) >
  //


  int      ix;
  int      iy;
  int      rxId;
  ifstream theFile;

  theFile.open( theFileName.c_str(), ios::in);

  int iread = 0;
  theFile >> ix >> iy >> rxId;
  while ( ! theFile.eof() ) {
    iread++;
    cout << ix <<  " " <<  iy << " " << rxId << endl;
    geomMap.insert( make_pair(rxId, make_pair( ix, iy) ) );
    theFile >> ix >> iy >> rxId;
  }

  cout << "FluxAnalyzerGen::readGeomFile:  read " << iread << " geometry lines."<< endl;

  theFile.close();
  return;
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3