[Bio] / RecoElemSBW / SbwFactory.cc Repository:
ViewVC logotype

View of /RecoElemSBW/SbwFactory.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (download) (as text) (annotate)
Tue Sep 23 20:22:29 2003 UTC (16 years, 2 months ago) by efrank
Branch: MAIN
CVS Tags: HEAD
Changes since 1.5: +14 -8 lines
SbwFactory:
	more fixes for writing SBML.

fixSbmlForCytoScape:

#
# libSbml-2.0.1 generates level 2 sbml.  Cytoscape only reads level 1.
# This sed script doesn't really convert 2 to 1.  It does enough hacks
# on the sbml to make it possible for cytoscape to read it.
#
# A better approach all around would be to generate a SIF file.
#

//***********************************************************************
// Reconstruction Readout Interface via SBW.
//
// Description:
//
//    Implementation comments only.  Generic comments go in .hh
//
//
// Author List:
//    Kaitlyn Hwang, Ed Frank, efrank@mcs.anl.gov
//
// History:
//    07 Juoy 03  efrank, hwang        First Version
// 
// BUGS:
//
//***********************************************************************

#include <stdio.h>
#include <string>
#include <vector>
#include <iterator>
#include <fstream>
#include <iostream>

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

#include "sbml/SBMLReader.h"
#include "sbml/SBMLTypes.h"
#include "sbml/SBMLWriter.h"
#include "SBW.h"

using namespace std;
using namespace SystemsBiologyWorkbench;


//************************************************************************
SbwFactory::SbwFactory( ) 
{}

//************************************************************************
RecoElem*
SbwFactory::getRecoElemFromSbmlFile( const std::string& filename ) {
//************************************************************************
  return getRecoElemFromSbmlFile( filename.c_str() );
}


//************************************************************************
RecoElem*
SbwFactory::getRecoElemById( unsigned int procId) {
//************************************************************************

  DataBlockWriter  args;
  DataBlockReader  result;

 
  SBW::connect();

  Module recoElemReadModule   = SBW::getModuleInstance("gov.mcs.anl.compbio");
  Service recoElemReadService = recoElemReadModule.findServiceByName("RecoElemRead");
  ServiceMethod RecoMethod    = recoElemReadService.getMethod("{}getTree(int)");

  // build arg list to call service:  pass it the node Id we want to get back
  // then actually do the call


  args << (int) procId;    // the int cast fixes an overload ambiguity
  result=RecoMethod.call( args );


  RecoElem* elem = sbwToRecoElem( result );

  recoElemReadModule.shutdown();

  SBW::disconnect();

  return elem;
}

//************************************************************************
void
SbwFactory::unpackRx( DataBlockReader& rIn, vector <Reaction>& rxns ) {
//************************************************************************

  int             id;
  int             dir;
  DataBlockReader inMetabL, outMetabL, enzymeL;
  DataBlockReader r;

  
  // The passed in DataBlockReader, rIn, begins with a List...shove it
  // into r and pull the rx info out of it.

  rIn >> r;
  
  r >> id;
  r >> dir;


  Reaction rxn( id, (Reaction::Direction) dir );
  
  vector<string>& inMetab   = rxn.getInputMetaboliteNames();
  vector<int>&    inMetabN  = rxn.getInputStoich();
  vector<string>& outMetab  = rxn.getOutputMetaboliteNames();
  vector<int>&    outMetabN = rxn.getOutputStoich();
  vector<string>& enzymes   = rxn.getEnzymeNames();
  int             nIn;
  int             nOut;
  int             nEnz;
  string          name;
  int             stoich;


  r >> inMetabL;
  inMetabL >> nIn;
  for ( int i=0; i<nIn; i++ ) {
    inMetabL >> name;
    inMetabL >> stoich;


    inMetab.push_back( name);
    inMetabN.push_back( stoich );
    
  
  }
 

  r >> outMetabL;
  outMetabL >> nOut;
  for ( int i=0; i<nOut; i++ ) {
    outMetabL >> name;
    outMetabL >> stoich;
    outMetab.push_back( name);
    outMetabN.push_back( stoich );
  }

  r >> enzymeL;
  enzymeL >> nEnz;
  for ( int i=0; i<nEnz; i++ ) {
    enzymeL >> name;
    enzymes.push_back( name );
  }

  

  rxns.push_back( rxn );
  

  return;
}

//************************************************************************
RecoElem*
SbwFactory::sbwToRecoElem( DataBlockReader& r ) {
//************************************************************************

  string          n;
  int             id;

  r >> n;
  r >> id;
  RecoElem* elem = new RecoElem( n, id );

  vector<Reaction>& rxns      = elem->getReactions();
  vector<RecoElem>& subElems  = elem->getSubElements();

  int nrx;
  r >> nrx;
  for( int i=0; i<nrx; i++ ) {
    unpackRx( r, rxns);
   
  }


  int nKid;
  r >> nKid;
  RecoElem* aKid;
  for ( int i=0; i<nKid; i++ ) {
    aKid = sbwToRecoElem( r );
    subElems.push_back( *aKid );
  }

 
  return elem;
}

//************************************************************************
RecoElem*
SbwFactory::getRecoElemFromSbmlFile( const char* filename ) {
//************************************************************************

 
  SBMLDocument_t *d;
  Model_t *m;

  d     = readSBML(filename);

  m = d->model;
  if (0==m) {
    cout << "SbmlFactory:: Did not read model " << filename << endl;
    return 0;
  }

  RecoElem* elem = SBMLtoRecoElem (m); 

  return elem;
}


//**********************************************************************
RecoElem*
SbwFactory::SBMLtoRecoElem(Model_t* m) {
//************************************************************************

  string name;
  int id, numreactant, numproduct, modNum;

  Reaction_t* r;  
  Species_t* s;
  ListOf_t* rlist = m->reaction;
 
  SpeciesReference_t* sr;
  

  //not sure what to put as an id, 0 for now

  if (0==m->name)   {
    name=m->id;
  }  else  {
    name=m->name;
  }
  cout << "model id/name is " << m->id << "/" << m->name << endl;  

  id=0;

  RecoElem* elem= new RecoElem( name, id);

  vector<Reaction>& rxns      = elem->getReactions();

  //put in reactions
  int reactionsize= ListOf_getNumItems(rlist);
  cout << "in SBMLFactory: number of reactions " << reactionsize << endl;
  

  for (int j=0; j< reactionsize ; j++)
    {
      
      //  cout << "Reaction : " << j <<endl;
      r=Model_getReaction(m,j);
  
      cout << "*****************" << endl;
      inreaction(m, r,rxns);

      cout << endl;
    }

  return elem;
  


}
//************************************************************************
void
SbwFactory::inreaction(Model_t* theModel, Reaction_t* r, vector<Reaction>& rxns) {
//************************************************************************

  int id;
  int dir;

  int numreactant=Reaction_getNumReactants(r);
  int numproduct=Reaction_getNumProducts(r);
  SpeciesReference_t* sr, *pr;
  Species_t*           species;
   
  //id not sure what to put here
  //direction is reversible  or not. false, true

  Reaction::Direction reDir;
  id=0;
  dir=r->reversible;
  Reaction rxn( id, Reaction::forward );
  rxn.setIdString( Reaction_getId( r ) );
  if ( dir ) {
    rxn.setDirection( Reaction::both );
  }

  vector<string>& reactant = rxn.getInputMetaboliteNames();
  vector<string>& product  = rxn.getOutputMetaboliteNames();
  vector<string>& enzymes  = rxn.getEnzymeNames();
  

  vector<int>& inMetabN = rxn.getInputStoich();
  vector<int>& outMetabN = rxn.getOutputStoich();
  
  
  //add reactants
  //add products
  
 

  //figuring out the common denominator
  int rdem = 1;
  int pdem = 1;
  int demstoich;

  for (int j=0; j< numreactant ; j++ ) {
      sr=Reaction_getReactant(r,j);
      rdem= rdem * sr->denominator;
  }

  
  for (int m=0; m<numproduct; m++) {
      pr=Reaction_getProduct(r,m);
      pdem= pdem * pr->denominator;
  }

  demstoich = rdem * pdem;
  //cout << "common denominator: " << demstoich;

   for (int k=0; k < numreactant ; k++ ) {
       sr=Reaction_getReactant(r,k);
       //cout << "test reactant: " << sr->species << endl;

       species = Model_getSpeciesById(theModel, sr->species);
       reactant.push_back ( Species_getName( species) );
       inMetabN.push_back( int( demstoich * sr->stoichiometry/sr->denominator));
   }


   for (int l=0; l<numproduct; l++) {

        pr=Reaction_getProduct(r,l);
	//cout << "product: " << pr->species  << endl;
	 
	species = Model_getSpeciesById(theModel, sr->species);
	reactant.push_back ( Species_getName( species) );

	//product.push_back(pr->species);
	outMetabN.push_back(int( demstoich * pr->stoichiometry/ pr->denominator ));
   }
   
   getReactionModifier(r, enzymes);

   rxns.push_back(rxn);
  
  return;
}

//************************************************************************
void
SbwFactory::getReactionModifier(Reaction_t* r, vector<string>& enzymeList ) {
//************************************************************************

  ModifierSpeciesReference_t* rMod;
  int modNum;


  modNum = Reaction_getNumModifiers(r);
  cout << modNum << " modifiers for name/id " << r->name << "/" << r->id << endl;

  for (int k =0; k < modNum ; k++)
    {
      rMod =Reaction_getModifier(r, k);
      cout << "the modifier name is " <<  ModifierSpeciesReference_getSpecies(rMod) << endl;
      enzymeList.push_back(ModifierSpeciesReference_getSpecies(rMod) );
    }
}

//************************************************************************
void
SbwFactory::writeSbmlFile( RecoElem& re, const char* filename) {
//************************************************************************
  SBMLDocument_t * document;
  
  document = SBMLDocument_create (); 
  document->model = Model_createWithName( re.getName().c_str() );
  
  buildSbmlModelFromRecoElem( document->model, re);
  writeSBML(document, filename);

  return;
}

//************************************************************************
void
SbwFactory::buildSbmlModelFromRecoElem( Model_t* theModel, 
					RecoElem& re) {
//************************************************************************

  std::set<Reaction> rxns;
  std::set<string>   metabolites;
  std::set<string>   enzymeNames;
  Reaction_t*        theSbmlRx;

  re.gatherReactants( rxns, metabolites, enzymeNames );

  // we want one SBML species definition per reactant even if that
  // reactant occurs in multiple RX.  So build a map from metabolite
  // and enzyme names to SBML Species.

  map<string, Species_t*> speciesMap;
  buildSpeciesMap( speciesMap, metabolites );
  buildSpeciesMap( speciesMap, enzymeNames );

  std::map<string, Species_t*>::iterator iSpecies    = speciesMap.begin();
  std::map<string, Species_t*>::iterator iSpeciesEnd = speciesMap.end();
  for( ; iSpecies != iSpeciesEnd; iSpecies++) {
    Model_addSpecies(theModel, iSpecies->second );
  }

  std::set<Reaction>::iterator  iRxn    = rxns.begin();
  std::set<Reaction>::iterator  iRxnEnd = rxns.end();
  for( ; iRxn != iRxnEnd; iRxn++ ) {
    cout << "[" << iRxn->getIdString() << "] " << iRxn->asString() << endl;
    theSbmlRx = Model_createReaction( theModel );
    Reaction_setId( theSbmlRx, iRxn->getIdString() );
    Reaction_setName( theSbmlRx, iRxn->getIdString() );
    if (iRxn->getDirection() == Reaction::both ) {
      Reaction_setReversible( theSbmlRx, 1 );
    } else if (iRxn->getDirection() == Reaction::backward ) {
      cout << "WARNING: backwards RX not handled correctly" << endl;
      Reaction_setReversible( theSbmlRx, 1 );
    } else {
      Reaction_setReversible( theSbmlRx, 0 );
    }
    
    addInOuts( iRxn->getInputMetaboliteNames(), 
	       iRxn->getInputStoich(), 
	       theModel, speciesMap,
	       &Model_createReactant );
    addInOuts( iRxn->getOutputMetaboliteNames(), 
	       iRxn->getOutputStoich(), 
	       theModel, speciesMap,
	       &Model_createProduct );
    addModifiers( iRxn->getEnzymeNames(), theModel, speciesMap );
  }

  return;
}

//************************************************************************
void
SbwFactory::addInOuts( const vector<string>&    names,
		       const vector<int>&       stoich,
		       Model_t*                 model,
		       map<string, Species_t*>& speciesMap,
		       SpeciesReference_t* (*createFn)(Model_t*) ) {
//************************************************************************
  Species_t* theSpecies;
  SpeciesReference_t* theSpeciesRef;

  int n=names.size();
  for( int i=0; i<n; i++) {
    theSpecies = speciesMap[ names[i] ];
    //    theSpeciesRef = Model_createReactant( model );
    theSpeciesRef = createFn( model );
    SpeciesReference_setSpecies(theSpeciesRef, Species_getName(theSpecies));
    SpeciesReference_setStoichiometry(theSpeciesRef, stoich[i] );
    SpeciesReference_setDenominator(theSpeciesRef, 1);
  }
}


//************************************************************************
void
SbwFactory::addModifiers( const vector<string>&    names,
			  Model_t*                 model,
			  map<string, Species_t*>& speciesMap) {
//************************************************************************
  Species_t* theSpecies;
  ModifierSpeciesReference_t* theSpeciesRef;

  int n=names.size();
  for( int i=0; i<n; i++) {
    theSpecies = speciesMap[ names[i] ];
    theSpeciesRef = Model_createModifier(model);
    ModifierSpeciesReference_setSpecies(theSpeciesRef, 
					Species_getName(theSpecies));
  }
}

//************************************************************************
void
SbwFactory::buildSpeciesMap( map<string, Species_t*>& speciesMap,
			     set<string>& names ) {
//************************************************************************

  Species_t*                   theSpecies;
  std::set<string>::iterator   iName         = names.begin();
  std::set<string>::iterator   iNameEnd      = names.end();
  
  for ( ; iName != iNameEnd; iName++ ) {
    theSpecies = Species_create();
    Species_initDefaults( theSpecies );
    Species_setInitialAmount(theSpecies,  0.0 );
    Species_setName( theSpecies, iName->c_str() );
    //Species_setId( theSpecies, iName->c_str() );     //$$$ WRONG
    speciesMap.insert(  make_pair( *iName, theSpecies) );
  }
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3