rk.cpp

Go to the documentation of this file.

#include <sstream>
#include <assert.h>
#include "rk.h"

RK::RK() {

}
RK::RK(const RK &rk) {

        this->knock_reacs=rk.knock_reacs;
        this->key=rk.key;
        this->metab_network=rk.metab_network;

        this->inactive_reacs=rk.inactive_reacs;
        this->uncollapsed_ineff_efm_indices.assign(rk.uncollapsed_ineff_efm_indices.begin(),rk.uncollapsed_ineff_efm_indices.end());
        this->uncollapsed_eff_efm_indices.assign(rk.uncollapsed_eff_efm_indices.begin(),rk.uncollapsed_eff_efm_indices.end());
}
RK::RK(map<int,string> knock_reacs, MetabolicNetwork *metab_network) {


        stringstream os;
        map<int,string>::iterator pi;

        //find key for the RK instance
        //sort pi->first
        for (pi=knock_reacs.begin();pi!=knock_reacs.end();pi++) {
                if (pi!=knock_reacs.begin())
                        os<<".";
                os << pi->first;
                //this->knock_reacs_nos.push_back(pi->first);
        }
        this->key= os.str();
        
        this->metab_network=metab_network;


        this->uncollapsed_ineff_efm_indices.assign(metab_network->get_ineff_efm_indices().begin(),metab_network->get_ineff_efm_indices().end());
        this->uncollapsed_eff_efm_indices.assign(metab_network->get_eff_efm_indices().begin(),metab_network->get_eff_efm_indices().end());

        //add each reaction to the current instance from knock_reacs
        add_reacs(knock_reacs);

}

int RK::add_reacs(map<int,string> reacs) {
        //read in local variables the EFM matrix, and indices of inefficient modes
        vector<string> reacs_names=this->metab_network->get_reacs_names();
        
        //update remain_ineff_modes_cnt;
        //in one scan find modes that have 0 flux for reactions in this->knock_reacs but non-zero for reaction reac
        //find how many modes will be collapsed with the deletion of reaction reac

        map<int,string>::iterator pi;   

        vector<long> remain_efm_indices;
        
        for (long i=0;i< this->uncollapsed_ineff_efm_indices.size();i++)  {

                int collapsed=1;
                for (pi=reacs.begin();pi!=reacs.end();pi++) {

                        int is_zero;
                        this->metab_network->is_efm_elem_zero(this->uncollapsed_ineff_efm_indices[i],pi->first, is_zero); 
                        if (!is_zero)  {
                                collapsed=0;
                                break;
                        }
                }
                if (collapsed)
                        remain_efm_indices.push_back(this->uncollapsed_ineff_efm_indices[i]);                   
        }
        this->uncollapsed_ineff_efm_indices.assign(remain_efm_indices.begin(),remain_efm_indices.end());
        remain_efm_indices.resize(0);

        for (long i=0;i< this->uncollapsed_eff_efm_indices.size();i++)  {

                int collapsed=1;
                for (pi=reacs.begin();pi!=reacs.end();pi++) {

                        int is_zero;
                        this->metab_network->is_efm_elem_zero(this->uncollapsed_eff_efm_indices[i],pi->first, is_zero); 
                        if (!is_zero)  {
                                collapsed=0;
                                break;
                        }
                }
                if (collapsed)
                        remain_efm_indices.push_back(this->uncollapsed_eff_efm_indices[i]);                     
        }
        this->uncollapsed_eff_efm_indices.assign(remain_efm_indices.begin(),remain_efm_indices.end());
        remain_efm_indices.resize(0);


        //now find what fraction of modes indexed by remain_ineff_efm_indices are collapsed if  'reac' is deleted

        
        for (pi=reacs.begin();pi!=reacs.end();pi++) {
                this->knock_reacs[pi->first]=reacs_names[pi->first];
        }

        //recompute key from scratch
        stringstream os;

        //find key for the RK instance
        //sort pi->first
        for (pi=knock_reacs.begin();pi!=knock_reacs.end();pi++) {
                if (pi!=knock_reacs.begin())
                        os<<".";
                os << pi->first;
        
        }


        this->key= os.str();
        

        return 0;
}


int RK::add_reac_direct(int reac,double &rank) {

        //read in local variables the EFM matrix, and indices of inefficient modes
        vector<string> reacs_names=this->metab_network->get_reacs_names();

        //in one scan find modes that have 0 flux for reactions in this->knock_reacs but non-zero for reaction reac
        //find how many modes will be collapsed with the deletion of reaction reac
        vector<long> remain_efm_indices;

        //now find what fraction of modes indexed by remain_ineff_efm_indices are collapsed if  'reac' is deleted
        long num_remain_ineff_efm_indices, num_ineff_collapsed=0;
        for (long i=0;i<this->uncollapsed_ineff_efm_indices.size();i++)  {
                int is_zero;
                this->metab_network->is_efm_elem_zero(this->uncollapsed_ineff_efm_indices[i],reac, is_zero); 
                if (is_zero)  
                        remain_efm_indices.push_back(this->uncollapsed_ineff_efm_indices[i]);
        }
        num_ineff_collapsed=uncollapsed_ineff_efm_indices.size()-remain_efm_indices.size();
        num_remain_ineff_efm_indices=remain_efm_indices.size();
        this->uncollapsed_ineff_efm_indices.assign(remain_efm_indices.begin(),remain_efm_indices.end());


        remain_efm_indices.resize(0);
        //now find what fraction of modes indexed by remain_ineff_efm_indices are collapsed if  'reac' is deleted
        long num_remain_eff_efm_indices,num_eff_collapsed=0;
        for (long i=0;i<this->uncollapsed_eff_efm_indices.size();i++)    {
                int is_zero;
                this->metab_network->is_efm_elem_zero(this->uncollapsed_eff_efm_indices[i],reac, is_zero); 
                if (is_zero)  
                        remain_efm_indices.push_back(this->uncollapsed_eff_efm_indices[i]);
        }
        num_eff_collapsed=uncollapsed_eff_efm_indices.size()-remain_efm_indices.size();
        num_remain_eff_efm_indices=remain_efm_indices.size();
        this->uncollapsed_eff_efm_indices.assign(remain_efm_indices.begin(),remain_efm_indices.end());

        if (num_eff_collapsed==0) {
                rank=-1.0;
        } else {

                if (this->metab_network->get_options().get_network_capacity().compare("high")==0) {

                        double frac_ineff_collapsed=num_ineff_collapsed*1.0/num_remain_ineff_efm_indices;
                        double frac_eff_collapsed=num_eff_collapsed*1.0/num_remain_eff_efm_indices;
                        rank=frac_ineff_collapsed*(1-frac_eff_collapsed);  //or frac_ineff_collapsed*frac_eff_collapsed

                } else if (this->metab_network->get_options().get_network_capacity().compare("low")==0) {

                        int not_all_eff_collapsed=(num_remain_eff_efm_indices!=num_eff_collapsed);
                        rank=((num_ineff_collapsed+num_eff_collapsed)*1.0/(num_remain_eff_efm_indices+num_remain_ineff_efm_indices))*not_all_eff_collapsed;  //or frac_ineff_collapsed*frac_eff_collapsed

                }


        }

        this->knock_reacs[reac]=reacs_names[reac];

        //recompute key from scratch
        stringstream os;

        map<int,string>::iterator pi;

        //find key for the RK instance
        //sort pi->first
        for (pi=knock_reacs.begin();pi!=knock_reacs.end();pi++) {
                if (pi!=knock_reacs.begin())
                        os<<".";
                os << pi->first;
        }

        this->key= os.str();

        return 0;
}

int RK::add_reac(int reac,long &rank) {

        //read in local variables the EFM matrix, and indices of inefficient modes
        vector<string> reacs_names=this->metab_network->get_reacs_names();
        
        //update remain_ineff_modes_cnt;
        //in one scan find modes that have 0 flux for reactions in this->knock_reacs but non-zero for reaction reac
        //find how many modes will be collapsed with the deletion of reaction reac

        
        vector<long> remain_efm_indices;

        for (long i=0;i< uncollapsed_ineff_efm_indices.size();i++)  {
                map<int,string>::iterator pi;
                int collapsed=1;

                int is_zero;
                this->metab_network->is_efm_elem_zero(uncollapsed_ineff_efm_indices[i],reac, is_zero); 
                if (is_zero)  
                        remain_efm_indices.push_back(uncollapsed_ineff_efm_indices[i]);                 
        }
        long num_collapsed;
        num_collapsed=uncollapsed_ineff_efm_indices.size()-remain_efm_indices.size();
        uncollapsed_ineff_efm_indices.assign(remain_efm_indices.begin(),remain_efm_indices.end());

        remain_efm_indices.resize(0);
        for (long i=0;i< uncollapsed_eff_efm_indices.size();i++)  {
                map<int,string>::iterator pi;
                int collapsed=1;

                int is_zero;
                this->metab_network->is_efm_elem_zero(uncollapsed_eff_efm_indices[i],reac, is_zero); 
                if (is_zero)  
                        remain_efm_indices.push_back(uncollapsed_eff_efm_indices[i]);                   
        }
        
        uncollapsed_eff_efm_indices.assign(remain_efm_indices.begin(),remain_efm_indices.end());

        //now find what fraction of modes indexed by remain_ineff_efm_indices are collapsed if  'reac' is deleted


        rank=num_collapsed; //*1.0/remain_ineff_efm_indices.size();

        this->knock_reacs[reac]=reacs_names[reac];

        //recompute key from scratch
        stringstream os;

        map<int,string>::iterator pi;

        //find key for the RK instance
        //sort pi->first
        for (pi=knock_reacs.begin();pi!=knock_reacs.end();pi++) {
                if (pi!=knock_reacs.begin())
                        os<<".";
                os << pi->first;
        }

        this->key= os.str();

        return 0;
}


//find the new rank value if reac is added
int RK::test_reac(int reac,long &rank) {

        //read in local variables the EFM matrix, and indices of inefficient modes

        vector<long>& ineff_efm_indices=this->metab_network->get_ineff_efm_indices();
        //update remain_ineff_modes_cnt;
        //in one scan find modes that have 0 flux for reactions in this->knock_reacs but non-zero for reaction reac
        //find how many modes will be collapsed with the deletion of reaction reac

        //store in remain_ineff_efm_indices the indices of those modes which are collapsed when reactions from reac_knocks are deleted
        vector<long> remain_ineff_efm_indices_lc;
        for (long i=0;i<ineff_efm_indices.size();i++)  {
                map<int,string>::iterator pi;
                int collapsed=1;
                for (pi=knock_reacs.begin();pi!=knock_reacs.end();pi++)  {
                        int is_zero;
                        this->metab_network->is_efm_elem_zero(ineff_efm_indices[i],pi->first, is_zero); 
                        if (!is_zero)  {
                                collapsed=0;
                                break;
                        }
                }
                if (collapsed)
                        remain_ineff_efm_indices_lc.push_back(ineff_efm_indices[i]);                    
        }

        //now find what fraction of modes indexed by remain_ineff_efm_indices are collapsed if  'reac' is deleted
        long num_collapsed=0;
        for (long i=0;i<remain_ineff_efm_indices_lc.size();i++) {
                int is_zero;
                this->metab_network->is_efm_elem_zero(remain_ineff_efm_indices_lc[i],reac, is_zero); 
                if (!is_zero)  
                        num_collapsed++;
        }

        rank=num_collapsed; //*1.0/remain_ineff_efm_indices.size();

        
        return 0;
}

long RK::test_knock_reacs2(map<int,string> knock_reacs, vector<double> &min_yields, long &num_eff_efm_rem, long &num_ineff_efm_rem) {


        num_ineff_efm_rem=this->uncollapsed_ineff_efm_indices.size();
        num_eff_efm_rem=this->uncollapsed_eff_efm_indices.size();

        vector< vector<double> >& efm_yields=this->metab_network->get_efm_yields();
        //find max yield for all target reactions among modes in remain_ineff_efm_indices
        int target_cnt=0;
        map<string,double>& max_yields=this->metab_network->get_max_yields();
        for (map<string,double>::iterator it=max_yields.begin();it!=max_yields.end();it++,target_cnt++) {
                //it->second is the index of yield reaction in efms
                double curr_min_yield=it->second;
                vector<long>::iterator ip;
                for (ip=this->uncollapsed_eff_efm_indices.begin();ip!=this->uncollapsed_eff_efm_indices.end();ip++) {
                        //*ip index of mode in efms
                        if (efm_yields[target_cnt][*ip]<curr_min_yield) 
                                curr_min_yield=efm_yields[target_cnt][*ip];
                
                }
                min_yields.push_back(curr_min_yield);
        }

        vector<long>::iterator ip;
        vector<string>& reacs_names=this->metab_network->get_reacs_names();
        this->max_path_length=0;this->min_path_length=reacs_names.size();


        this->avg_path_length=0;

        for (ip=this->uncollapsed_eff_efm_indices.begin();ip!=this->uncollapsed_eff_efm_indices.end();ip++) {
                int curr_path_length=0;
                for (int i=0;i<reacs_names.size();i++) {
                        int is_zero;
                        this->metab_network->is_efm_elem_zero(*ip,i, is_zero); 
                        if (!is_zero)
                                curr_path_length++;
                }
                if (this->min_path_length>curr_path_length)
                        this->min_path_length=curr_path_length;
                if (this->max_path_length<curr_path_length)
                        this->max_path_length=curr_path_length;
                this->avg_path_length+=curr_path_length;
        }

        for (int i=0;i<reacs_names.size();i++) {
                //determine if i-th reaction is inactive, but not in knock_reacs
                map<int,string>::iterator mis=knock_reacs.find(i);
                if (mis!=knock_reacs.end())
                        continue;
                int all_zero=1;
                for (ip=this->uncollapsed_eff_efm_indices.begin();ip!=this->uncollapsed_eff_efm_indices.end();ip++) {
                        int is_zero;
                        this->metab_network->is_efm_elem_zero(*ip,i, is_zero); 
                        if (!is_zero)   {
                                all_zero=0;
                        }               
                }
                if (all_zero)
                        this->inactive_reacs.push_back(i);
        }
        assert(this->uncollapsed_eff_efm_indices.size()>0);
        
        this->avg_path_length/=this->uncollapsed_eff_efm_indices.size();
        return 0;
}

ostream& operator<<(ostream& out, const RK& rk) {
        RK temp_rk=rk;
        map<string,int>  transport_names=rk.metab_network->get_transport_names();

        for (map<int,string>::iterator it=temp_rk.get_knock_reacs().begin();it!=temp_rk.get_knock_reacs().end();it++) {
                //check if it->second is transport reaction
                
                map<string,int>::iterator mpi=transport_names.find(it->second);
                if (mpi==transport_names.end()) {
                        out << it->second << " ";
                } else
                        out << "t"<<it->second << " ";
        }
        //find inactive reactions in all modes remain_eff_efm_indices modes and print those which don't appear in rk.get_knock_reacs()

        out << "\t";
        
        vector<string>& reacs_names=rk.metab_network->get_reacs_names();
        for (int i=0;i<rk.inactive_reacs.size();i++) {
                out << reacs_names[rk.inactive_reacs[i]]<<" ";
        }
        
        return out;

}
int RK::is_core_collapsed(int &cnt_not_collapsed) {
        //does knock_reacs collapse all modes in this->metab_network->get_core_efm_indices();   

        vector< vector<double> >& core_efms=this->metab_network->get_core_efms();
        map<int,string>::iterator it;
        cnt_not_collapsed=0;

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

        }

        return cnt_not_collapsed==0;
}

int find_common_efm(vector<RK> rks, int &num_shared_efms) {

        
        vector<long> common_efm=rks[0].get_uncollapsed_eff_efm_indices();       

        
        for (vector<RK>::iterator it=rks.begin()+1;it!=rks.end();it++)  {
                vector<long> remain_efm_ind=(*it).get_uncollapsed_eff_efm_indices();
                sort(remain_efm_ind.begin(),remain_efm_ind.end());
                sort(common_efm.begin(),common_efm.end());
                vector<long> vec2;
                std::set_intersection(common_efm.begin(),common_efm.end(),remain_efm_ind.begin(),remain_efm_ind.end(),std::back_inserter(vec2));
                common_efm=vec2;
        }
        
        num_shared_efms=common_efm.size();
        return 0;
}