metabolicnetwork.cpp

Go to the documentation of this file.

#include <iostream>
#include <fstream>   // file I/O
#include <sstream>
#include <stdlib.h>
#include <vector>
#include <algorithm>
#include <math.h>
#include <iomanip>
#include <limits>
#include "metabolicnetwork.h"
#include "util.h"
#include "mpi_utils.h"



MetabolicNetwork::MetabolicNetwork() {}
MetabolicNetwork::MetabolicNetwork(MetabolicNetwork &metab_network) {


        this->efm_yields=efm_yields;

        //init reacs_names
        this->reacs_names=metab_network.reacs_names;
        this->substrate_names=metab_network.substrate_names;
        //init target_names
        this->target_names=metab_network.target_names;
        this->yields=metab_network.yields;
        this->num_modes=metab_network.num_modes;
        this->num_reacs=metab_network.num_reacs;
        this->num_metabs=metab_network.num_metabs;


}

MetabolicNetwork::MetabolicNetwork(Options options):bit_efm()  {

        vector< vector<double> > stoich;
        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);
        t_UnsignedInt num_ems;

        this->options=options;

        this->options.load_efm_meta_file(stoich, this->reacs_names, num_ems);

        ifstream fp_in;  // file for reading
        ofstream fp_out;

        vector<string> init_knock_reacs=options.get_init_knock_reacs();
        vector<string> target_reacs=options.get_target_reacs();
        vector<string> substrate_reacs=options.get_substrate_reacs();
        vector<double> min_yields=options.get_min_yields();
        vector<double> min_core_yields=options.get_min_core_yields();
        vector<string> transport_reacs=options.get_transport_reacs();


        this->compressed_mode=1;
        bit_efm.load_efm_bit(options.get_efm_bit_fname(), num_ems, stoich);

        /*****************INIT DATA STRUCTURES FOR WHAT IS LOADED FROM configs.knock******************/
        //read substrate reactions
        for (vector<string>::iterator it=substrate_reacs.begin();it!=substrate_reacs.end();it++) {
                if ((*it).size()) {
                    vector<string>::iterator reac_ind=find(this->reacs_names.begin(),this->reacs_names.end(),*it);
                    this->substrate_names[*it]=reac_ind-reacs_names.begin();
                }               
        }


        map<string,int>::iterator pi;
        //if (proc_id==0) {
        //      cout << endl << "\t"<<util::get_time_str()<<": substrates: ";
        //      for (pi=substrate_names.begin();pi!=substrate_names.end();pi++)
        //              cout << " "<<pi->first;
        //      cout << endl;
        //}


        //read transport reactions
        for (vector<string>::iterator it=transport_reacs.begin();it!=transport_reacs.end();it++) {
                if ((*it).size()) {
                    vector<string>::iterator reac_ind=find(this->reacs_names.begin(),this->reacs_names.end(),*it);
                    this->transport_names[*it]=reac_ind-reacs_names.begin();
                }               
        }
        
        //if (proc_id==0) {
        //      cout << endl << "\t"<<util::get_time_str()<<": transport reactions: ";
        //      for (pi=transport_names.begin();pi!=transport_names.end();pi++)
        //              cout << " "<<pi->first;
        //      cout << endl;
        //}

        
        if (options.get_use_xchg_reacs().compare("yes")) {
                //remove transport reacs from knock_cand_reacs
                int i=0;
                for (vector<string>::iterator it=this->reacs_names.begin();it!=this->reacs_names.end();it++,i++) {
                        //is reaction *it in transport_names
                        map<string,int>::iterator mpi=this->transport_names.find(*it);
                        if (mpi==this->transport_names.end()) {
                                knock_cand_reacs.push_back(i);
                        }
                }
        } else {
                for (int i=0;i<this->reacs_names.size();i++)
                        knock_cand_reacs.push_back(i);
        }


        if (proc_id==0) {
                cout << endl << "\t"<<util::get_time_str()<<": knockout candidate reactions: ";
                for (int i=0;i<this->knock_cand_reacs.size();i++)
                        cout << " "<<this->reacs_names[this->knock_cand_reacs[i]];
                cout << endl;
        }
                
        
        //read target reactions
        
        for (vector<string>::iterator it=target_reacs.begin();it!=target_reacs.end();it++) {
                if ((*it).size()) {
                        vector<string>::iterator reac_ind;
                        reac_ind=find(reacs_names.begin(),reacs_names.end(),*it);
                        this->target_names[*it]=reac_ind-reacs_names.begin();
                }               
        }
        
        //read init knockout reactions (e.g. for anaerobic case OPM1, OPM2)

        for (vector<string>::iterator it=init_knock_reacs.begin();it!=init_knock_reacs.end();it++) {

                if ((*it).size()) {
                        vector<string>::iterator reac_ind;
                        reac_ind=find(reacs_names.begin(),reacs_names.end(),*it);
                        this->init_knocks_names[*it]=reac_ind-reacs_names.begin();
                }               
        }

        //if (proc_id==0) {
        //      cout << endl << "\t"<<util::get_time_str()<<": initial knockout reactions"; // << this->init_knocks_names.size()<<endl;
        //      map<string,int>::iterator t;
        //      for (t=this->init_knocks_names.begin();t!=this->init_knocks_names.end();t++)
        //              cout <<" ("<< t->first<<","<<t->second<<")";
        //      cout << endl;
        //}
                //read target min yield values

        map<string,int>::iterator p=this->target_names.begin();

        for (vector<double>::iterator dit=min_yields.begin();dit!=min_yields.end();dit++) {
                this->yields[p->first]=*dit; //yield;
                p++;            
        }


        map<string,double>::iterator q;

        if (proc_id==0) {
                cout << endl <<"\t"<<util::get_time_str()<<": target reactions: ";
                for (p=this->target_names.begin() ,q=this->yields.begin();p!=this->target_names.end();p++,q++ )
                        cout << " (" << p->first <<","<<p->second<<")|"<<q->second;             
        }
        /***********************************************************************/       


        init_yields();

        remove_zero_substrate_modes();

        if (proc_id==0)
                cout << endl << "\t"<<util::get_time_str()<<": # of EFMs after removing modes with zero-flux substrates "<< get_num_modes()<<endl;
        
        bit_efm.knock_efms(this->init_knocks_names, this->efm_yields);

        if (proc_id==0)
                cout << endl << "\t"<<util::get_time_str()<<": # of EFMs after initial knockout "<< get_num_modes();




        this->num_reacs=reacs_names.size();
        this->num_modes=get_num_modes();        

        

        if (proc_id==0) {
                cout << endl << "\t"<<util::get_time_str()<<": # of elementary modes: " << this->num_modes;
                cout << endl << "\t"<<util::get_time_str()<<": # of reactions: " << this->num_reacs;
        }
        
        //*******************************
        //read target min core EFM yield values

        init_max_yields();
        p=this->target_names.begin();
        if (proc_id==0)
                cout << endl <<"\t"<<util::get_time_str()<<": minimum production-efficient mode yields [0-1]: ";
        
        for (vector<double>::iterator dit=min_core_yields.begin();dit!=min_core_yields.end();dit++) {   
                if (proc_id==0)
                        cout << "("<< p->first<<","<<(*dit)*(this->max_yields[p->first])<<") ";
                this->core_min_yields.push_back((*dit)*(this->max_yields[p->first]));   //multiply core_yield with max_yield
                p++;            
        }

        if (proc_id==0)
                cout << endl;
        //end read target min core EFM yield values
        
        //separate core modes from efms
        extract_production_eff_efms();
        //find indices of efficient and inefficient modes
        extract_growth_eff_and_ineff_modes();
        
        
        if (proc_id==0) {
                cout << endl << "\t"<<util::get_time_str()<<": # of production-efficient modes: " << this->core_efm_indices.size();
                cout << endl << "\t"<<util::get_time_str()<<": # of growth-efficient modes: " << this->eff_efm_indices.size();
                cout << endl << "\t"<<util::get_time_str()<<": # of inefficient modes: " << this->ineff_efm_indices.size();
        }

}



int MetabolicNetwork::init_yields()  {
        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);


        if (proc_id==0) {
                fstream fin;
                fin.open(this->options.get_yields_fname().c_str(),ios::in);
                if( fin.is_open() ) {

                        cout << endl << "\t"<<util::get_time_str()<<": reading mode yields from file "<<this->options.get_yields_fname().c_str()<<".... ";
                        string str;
                                this->efm_yields.resize(this->target_names.size(),vector<double>(get_num_modes()));
                                for (long i=0;i<get_num_modes();i++)
                                        for (long j=0;j<this->target_names.size();j++) {
                                                fin >> str;
                                                if (str.compare("inf")==0)
                                                        this->efm_yields[j][i]=numeric_limits<double>::max();
                                                else if (str.compare("-inf")==0) 
                                                        this->efm_yields[j][i]=numeric_limits<double>::min();
                                                else {
                                                        istringstream ss(str);
                                                        ss>>this->efm_yields[j][i];
                                                }
                                        }
                }
                fin.close();
        }
        mpi_bcast_matrix( this->efm_yields);

        return 0;
}
int MetabolicNetwork::init_max_yields()  {

        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);

        map<string,int>::iterator p;
        for (p=this->target_names.begin();p!=this->target_names.end();p++) {
                this->max_yields[p->first]=0;
        }
        
        for (long j=0;j<this->efm_yields[0].size();j++) {
                int k=0;
                
                for (p=this->target_names.begin();p!=this->target_names.end();p++) {
        
                        if (this->efm_yields[k][j]>this->max_yields[p->first]) 
                                this->max_yields[p->first]=this->efm_yields[k][j];
                        k++;
                }
        
        }
        if (proc_id==0) {
                map<string,double>::iterator q,t;
                //cout << endl << "\t"<<util::get_time_str()<<": target reaction yields: ";
                //for (q=this->max_yields.begin(),t=this->yields.begin();q!=this->max_yields.end();q++,t++) {
                //      cout  << "("<< q->first << ","<<q->second<<","<<t->second<<")";
                //}
                cout << endl << "\t"<<util::get_time_str()<<": minimum growth-efficient mode yields [0-1]:    ";
                for (q=this->max_yields.begin(),t=this->yields.begin();q!=this->max_yields.end();q++,t++) {
                        cout  << "("<< q->first << ","<<q->second*t->second<<") ";
                }
        }

        
}

int MetabolicNetwork::get_num_reacs() { return this->reacs_names.size(); }

        
unsigned long MetabolicNetwork::get_num_modes() { 
                unsigned long val;
                val=this->bit_efm.get_num_modes();
                return val;
} 


int MetabolicNetwork::is_efm_elem_zero(unsigned long efm_index, int reac, int &is_zero) {


        this->bit_efm.is_efm_elem_zero(efm_index,reac,is_zero);
        return 0;

}


double MetabolicNetwork::get_efm_elem(unsigned long index, int reac) {
        double elem;
        vector<double> em;
        this->bit_efm.read_single_efm(index,em);
        elem=em[reac];
        return elem;
}


int MetabolicNetwork::remove_zero_substrate_modes() {

        //remove elementary modes having zero-flux for all substrate reactions
        vector<long> efm_zero_substrates;

        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);

        
        for (long j=0;j<get_num_modes();j++) {


                //now iterate through substrates
                map<string,int>::iterator subst_it;
                int all_subst_zero=1;
                for (subst_it=this->substrate_names.begin();subst_it!=this->substrate_names.end();subst_it++)  {

                        int is_zero;
                        is_efm_elem_zero(j, subst_it->second, is_zero); 

                        if (!is_zero) {
                                all_subst_zero=0;
                                break;
                        }
                }

                if (all_subst_zero==1)
                        efm_zero_substrates.push_back(j);
        }
        if (proc_id==0)
                cout << endl << "\t"<<util::get_time_str()<<": # of modes with zero-flux for substrate reaction: "<<efm_zero_substrates.size();

        this->bit_efm.remove_efm(efm_zero_substrates,this->efm_yields); 

        if (proc_id==0)
                cout << endl << "\t"<<util::get_time_str()<<": removed modes with 0 substrate flux";

}
int MetabolicNetwork::is_core_efms_collapsed(map<int,string> knock_reacs,int &cnt_not_collapsed) {

        map<int,string>::iterator it;
        cnt_not_collapsed=0;

        for (long i=0;i<this->core_efms[0].size();i++) {
                int not_collapsed=1;
                for (it=knock_reacs.begin();it!=knock_reacs.end();it++) {
                        if (this->core_efms[it->first][i]!=0) {
                                not_collapsed=0;
                                break;
                        }
                }
                if (not_collapsed)
                        cnt_not_collapsed++;

        }

        return cnt_not_collapsed==0;
}
int MetabolicNetwork::extract_production_eff_efms() {

        int nump=1,proc_id=0;
        util::get_proc_id_nump(proc_id,nump);

        for (long i=0;i<get_num_modes();i++) {
                //yield for efms[][j]-th mode is in array this->efm_yields[j][i]=
                int is_core_efm=1;
                for (int j=0;j<this->core_min_yields.size();j++)
                        if (this->efm_yields[j][i]<core_min_yields[j]) {
                                is_core_efm=0;
                                break;
                        }
                if (is_core_efm)
                        this->core_efm_indices.push_back(i);    
        }

        //extract core modes from this->efms and create this->core_efms
        this->core_efms.resize(this->num_reacs,vector<double>(this->core_efm_indices.size()));
        long k=0;
        for (long j=0;j<this->core_efm_indices.size();j++) {
                //found core mode
                for (int i=0;i<this->num_reacs;i++) {
                        double elem;
                        elem=get_efm_elem(this->core_efm_indices[j], i);
                        this->core_efms[i][k]=elem; 
                }                       
                k++;
        }

        bit_efm.remove_efm(this->core_efm_indices,this->efm_yields);

}
int MetabolicNetwork::extract_growth_eff_and_ineff_modes() {

        for (long j=0;j<get_num_modes();j++) {
                //explore yield of i-th EFM
                //double substrate_sum=0;
                map<string,int>::iterator p;
                int k=0;
                for (p=this->target_names.begin();p!=this->target_names.end();p++) {
                        if (this->efm_yields[k++][j]<(this->yields[p->first])*(this->max_yields[p->first]))
                                break; 
                }
                if (p!=this->target_names.end())
                        this->ineff_efm_indices.push_back(j);
                else
                        this->eff_efm_indices.push_back(j);                     
        }
        return 0;
}