Home > RAVEN > checkProduction.m

checkProduction

PURPOSE ^

checkProduction

SYNOPSIS ^

function [notProduced, notProducedNames, neededForProductionMat,minToConnect,model]=checkProduction(model,checkNeededForProduction,excretionFromCompartments,printDetails)

DESCRIPTION ^

 checkProduction
   Checks which metabolites that can be produced from a model using the
   specified constraints.

   model                       a model structure
   checkNeededForProduction    for each of the metabolites that could not
                               be produced, include an artificial 
                               production reaction and calculate which new
                               metabolites that could be produced as en
                               effect of this (opt, default false)
   excretionFromCompartments   cell array with compartment ids from which
                               metabolites can be excreted (opt, default
                               model.comps)
   printDetails                print details to the screen (opt, default
                               true)

   notProduced                 cell array with metabolites that could not
                               be produced
   notProducedNames            cell array with names and compartments for
                               metabolites that could not be produced
   neededForProductionMat      matrix where n x m is true if metabolite n
                               allows for production of metabolite m
   minToConnect                structure with the minimal number of
                               metabolites that need to be connected in 
                               order to be able to produce all other 
                               metabolites and which metabolites each of
                               them connects
   model                       updated model structure with excretion
                               reactions added

   The function is intended to be used to identify which metabolites must
   be connected in order to have a fully connected network. It does so by
   first identifying which metabolites could have a net production in the
   network. Then it calculates which other metabolites must be able to
   have net production in order to have production of all metabolites in
   the network. So, if a network contains the equations A[external]->B, 
   C->D, and D->E it will identify that production of C will connect 
   the metabolites D and E.

   Usage: [notProduced, notProducedNames,neededForProductionMat,minToConnect,model]=...
           checkProduction(model,checkNeededForProduction,...
           excretionFromCompartments,printDetails)

   Rasmus Agren, 2013-02-06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [notProduced, notProducedNames, neededForProductionMat,minToConnect,model]=checkProduction(model,checkNeededForProduction,excretionFromCompartments,printDetails)
0002 % checkProduction
0003 %   Checks which metabolites that can be produced from a model using the
0004 %   specified constraints.
0005 %
0006 %   model                       a model structure
0007 %   checkNeededForProduction    for each of the metabolites that could not
0008 %                               be produced, include an artificial
0009 %                               production reaction and calculate which new
0010 %                               metabolites that could be produced as en
0011 %                               effect of this (opt, default false)
0012 %   excretionFromCompartments   cell array with compartment ids from which
0013 %                               metabolites can be excreted (opt, default
0014 %                               model.comps)
0015 %   printDetails                print details to the screen (opt, default
0016 %                               true)
0017 %
0018 %   notProduced                 cell array with metabolites that could not
0019 %                               be produced
0020 %   notProducedNames            cell array with names and compartments for
0021 %                               metabolites that could not be produced
0022 %   neededForProductionMat      matrix where n x m is true if metabolite n
0023 %                               allows for production of metabolite m
0024 %   minToConnect                structure with the minimal number of
0025 %                               metabolites that need to be connected in
0026 %                               order to be able to produce all other
0027 %                               metabolites and which metabolites each of
0028 %                               them connects
0029 %   model                       updated model structure with excretion
0030 %                               reactions added
0031 %
0032 %   The function is intended to be used to identify which metabolites must
0033 %   be connected in order to have a fully connected network. It does so by
0034 %   first identifying which metabolites could have a net production in the
0035 %   network. Then it calculates which other metabolites must be able to
0036 %   have net production in order to have production of all metabolites in
0037 %   the network. So, if a network contains the equations A[external]->B,
0038 %   C->D, and D->E it will identify that production of C will connect
0039 %   the metabolites D and E.
0040 %
0041 %   Usage: [notProduced, notProducedNames,neededForProductionMat,minToConnect,model]=...
0042 %           checkProduction(model,checkNeededForProduction,...
0043 %           excretionFromCompartments,printDetails)
0044 %
0045 %   Rasmus Agren, 2013-02-06
0046 %
0047 
0048 if nargin<2
0049     checkNeededForProduction=false;
0050 end
0051 
0052 if nargin<3
0053     excretionFromCompartments=model.comps;
0054 end
0055 
0056 if nargin<4
0057     printDetails=true;
0058 end
0059 
0060 %Add an exchange reaction for each metabolite in the allowed compartments
0061 %and see if it can carry a flux
0062 allowedMetIds=ismember(model.comps(model.metComps),excretionFromCompartments);
0063 allowedMetIndexes=find(allowedMetIds);
0064 [model addedRxns]=addExchangeRxns(model,'out',allowedMetIndexes);
0065 
0066 canProduce=haveFlux(model,10^-5,addedRxns);
0067 
0068 notProduced=find(~canProduce);
0069 minToConnect={};
0070 if checkNeededForProduction==true
0071     %For each of the metabolites that couldn't be produced allow uptake and check
0072     %which of the other metabolites that couldn't be produced that can be
0073     %produced
0074     neededForProductionMat=false(numel(notProduced));
0075     for i=1:numel(notProduced)
0076         %Add uptake for this metabolite
0077         if i>1
0078             %Reset last iteration
0079             model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i-1))=model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i-1))*-1;
0080         end
0081         %Change the production reaction to an uptake reaction
0082         model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i))=model.S(:,numel(model.rxns)-numel(addedRxns)+notProduced(i))*-1;
0083         
0084         %Test which of the metabolites that couldn't be produced that can
0085         %be produced now
0086         neededForProductionMat(i,:)=haveFlux(model,10^-5,addedRxns(notProduced));
0087     end
0088     %Calculate the smallest number of metabolites that must be connected to
0089     %make everything connected and return their names
0090     
0091     %The algorithm is relatively straight forward. It finds the metabolite
0092     %that connects the most unconnected metabolites (iteratively), adds it
0093     %and removes the now connected metabolites until all are connected.
0094     %This is not guaranteed to find the global minimum
0095     neededForProdTemp=neededForProductionMat;
0096     while 1==1
0097         %Get all metabolites a metabolite makes connected
0098         totalConnected=false(size(neededForProdTemp));
0099         for i=1:numel(notProduced)
0100             totalConnected(i,:)=neededForProdTemp(i,:);
0101             
0102             lastIter=0;
0103             while 1==1
0104                 [crap a crap]=find(neededForProdTemp(totalConnected(i,:),:));
0105                 totalConnected(i,a)=true;
0106                 if numel(a)==lastIter
0107                    break; %No more connections were possible
0108                 else
0109                     lastIter=numel(a); 
0110                 end
0111             end
0112         end
0113         [connections mostConnected]=max(sum(totalConnected,2));
0114         
0115         if connections>0            
0116             %Add the most connected metabolite to the list and remove all
0117             %metabolites that it's connected to
0118             metID=allowedMetIndexes(notProduced(mostConnected));
0119             entry=[model.metNames{metID},'[',model.comps{model.metComps(metID)},'] (connects ' num2str(connections) ' metabolites)'];
0120             minToConnect=[minToConnect;entry];
0121             neededForProdTemp(totalConnected(mostConnected,:),:)=false;
0122             neededForProdTemp(:,totalConnected(mostConnected,:))=false;
0123         else
0124             break;
0125         end
0126     end
0127 else
0128     neededForProductionMat=[];
0129 end
0130 
0131 notProducedNames=strcat(model.metNames(allowedMetIndexes(notProduced)),'[',model.comps(model.metComps(allowedMetIndexes(notProduced))),']');
0132 
0133 if printDetails==true
0134     fprintf('The following metabolites could not be produced:\n');
0135     [notProducedNamesTemp,perm]=sort(notProducedNames);
0136     
0137     if checkNeededForProduction==true
0138         neededForProdTemp=neededForProductionMat(:,perm);
0139         neededForProdTemp=neededForProdTemp(perm,:);
0140         fprintf('\tIf the production of a metabolite is dependent on some other metabolites then those are printed under the name\n\n');
0141     end
0142     for i=1:numel(notProducedNamesTemp)
0143         fprintf([notProducedNamesTemp{i} '\n']);
0144         neededForProdTemp(i,i)=false; %Not neat to do this here. Prevent printing itself
0145         if checkNeededForProduction==true
0146             enablesProduction=find(neededForProdTemp(:,i));
0147             if any(enablesProduction)
0148                 for j=1:numel(enablesProduction)
0149                     fprintf(['\t' notProducedNamesTemp{enablesProduction(j)} '\n']);
0150                 end
0151             end
0152         end
0153     end
0154 end
0155 end

Generated on Mon 06-Jan-2014 14:58:12 by m2html © 2005