Merge pull request #505 from SysBioChalmers/develop

RAVEN 2.8.4

Author: edkerk

core/addExchangeRxns.m

@@ -89,4 +89,7 @@
 if isfield(model,'rxnConfidenceScores')
     model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(J),1)];
 end
+if isfield(model,'rxnDeltaG')
+    model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(J),1)];
+end
 end

core/addMets.m

@@ -29,6 +29,9 @@
 %               metMiriams     cell array with MIRIAM structures (opt,
 %                              default [])
 %               metCharges     metabolite charge (opt, default NaN)
+%               metDeltaG      Gibbs free energy of formation at
+%                              biochemical standard condition in kJ/mole
+%                              (opt, default NaN)
 %               metNotes       cell array with metabolite notes as strings
 %                              (opt, default '')
 %   copyInfo    when adding metabolites to a compartment where it previously
@@ -261,6 +264,27 @@
     end
 end
 
+if isfield(metsToAdd,'metDeltaG')
+    if numel(metsToAdd.metDeltaG)~=nMets
+        EM='metsToAdd.metDeltaG must have the same number of elements as metsToAdd.mets';
+        dispEM(EM);
+    end
+    if ~isnumeric(metsToAdd.metDeltaG)
+        EM='metsToAdd.metDeltaG must be of type "double"';
+        dispEM(EM);
+    end
+    if ~isfield(newModel,'metDeltaG')
+        newModel.metDeltaG=NaN(numel(largeFiller),1);
+    end
+    newModel.metDeltaG=[newModel.metDeltaG;metsToAdd.metDeltaG(:)];
+else
+    %Add default
+    if isfield(newModel,'metDeltaG')
+        newModel.metDeltaG=[newModel.metDeltaG;NaN(numel(filler),1)];
+    end
+end
+
+
 if isfield(metsToAdd,'metNotes')
     metsToAdd.metNotes=convertCharArray(metsToAdd.metNotes);
     if numel(metsToAdd.metNotes)==1 && numel(metsToAdd.mets)>1
@@ -345,6 +369,9 @@
         if isfield(newModel,'metCharges')
             newModel.metCharges(I(i))=newModel.metCharges(J(i));
         end
+        if isfield(newModel,'metDeltaG')
+            newModel.metDeltaG(I(i))=newModel.metDeltaG(J(i));
+        end
     end
 end
 end

core/addRxns.m

@@ -52,6 +52,8 @@
 %                               model.comps) (opt, default {})
 %            rxnNotes           cell array with reaction notes (opt,
 %                               default '')
+%            rxnDeltaG          Gibbs free energy at biochemical standard
+%                               condition in kJ/mole (opt, default NaN)
 %            rxnReferences      cell array with reaction references (opt,
 %                               default '')
 %            rxnConfidenceScores   vector with reaction confidence scores
@@ -504,8 +506,6 @@
     %Fill with standard if it doesn't exist
     if ~isfield(newModel,'rxnConfidenceScores')
         newModel.rxnConfidenceScores=NaN(nOldRxns,1);
-        EM='Adding reactions with confidence scores without such information. All existing reactions will have confidence scores as NaNs';
-        dispEM(EM,false);
     end
     newModel.rxnConfidenceScores=[newModel.rxnConfidenceScores;rxnsToAdd.rxnConfidenceScores(:)];
 else
@@ -515,6 +515,23 @@
     end
 end
 
+if isfield(rxnsToAdd,'rxnDeltaG')
+    if numel(rxnsToAdd.rxnDeltaG)~=nRxns
+        EM='rxnsToAdd.rxnDeltaG must have the same number of elements as rxnsToAdd.rxns';
+        dispEM(EM);
+    end
+    %Fill with standard if it doesn't exist
+    if ~isfield(newModel,'rxnDeltaG')
+        newModel.rxnDeltaG=NaN(nOldRxns,1);
+    end
+    newModel.rxnDeltaG=[newModel.rxnDeltaG;rxnsToAdd.rxnDeltaG(:)];
+else
+    %Fill with standard if it doesn't exist
+    if isfield(newModel,'rxnDeltaG')
+        newModel.rxnDeltaG=[newModel.rxnDeltaG;NaN(nRxns,1)];
+    end
+end
+
 
 %***Start parsing the equations and adding the info to the S matrix The
 %mets are matched to model.mets

core/addRxnsGenesMets.m

@@ -116,21 +116,24 @@
     if isfield(sourceModel,'metFormulas')
         metsToAdd.metFormulas=sourceModel.metFormulas(metIdx);
     end
-    if isfield(sourceModel,'metCharge')
-        metsToAdd.metCharge=sourceModel.metCharge(metIdx);
+    if isfield(sourceModel,'metCharges')
+        metsToAdd.metCharges=sourceModel.metCharges(metIdx);
     end
     if isfield(sourceModel,'metMiriams')
         metsToAdd.metMiriams=sourceModel.metMiriams(metIdx);
     end
-    if isfield(sourceModel,'metFormulas')
-        metsToAdd.metFormulas=sourceModel.metFormulas(metIdx);
+    if isfield(sourceModel,'metNotes')
+        metsToAdd.metNotes=sourceModel.metNotes(metIdx);
     end
     if isfield(sourceModel,'inchis')
         metsToAdd.inchis=sourceModel.inchis(metIdx);
     end
     if isfield(sourceModel,'metSmiles')
         metsToAdd.metSmiles=sourceModel.metSmiles(metIdx);
     end
+    if isfield(sourceModel,'metDeltaG')
+        metsToAdd.metDeltaG=sourceModel.metDeltaG(metIdx);
+    end
 
     metsToAdd.compartments=strtrim(cellstr(num2str(sourceModel.metComps(metIdx)))); % Convert from compartment string to compartment number
     [~,idx]=ismember(metsToAdd.compartments,strsplit(num2str(1:length(sourceModel.comps)))); % Match compartment number to compartment abbreviation
@@ -164,6 +167,9 @@
     dispEM(EM, true);
 end
 rxnToAdd.rxnConfidenceScores(:)=confidence;
+if isfield(sourceModel,'rxnDeltaG')
+    rxnToAdd.rxnDeltaG=sourceModel.rxnDeltaG(rxnIdx);
+end
 if isfield(sourceModel,'subSystems')
     rxnToAdd.subSystems=sourceModel.subSystems(rxnIdx);
 end

core/addTransport.m

@@ -167,6 +167,9 @@
     if isfield(model,'rxnConfidenceScores')
         model.rxnConfidenceScores=[model.rxnConfidenceScores;ones(nRxns,1)];
     end
+    if isfield(model,'rxnDeltaG')
+        model.rxnDeltaG=[model.rxnDeltaG;zeros(nRxns,1)];
+    end
     addedRxns = [addedRxns; addedRxnsID];
 end
 end

core/changeRxns.m

@@ -129,6 +129,9 @@
 if isfield(model,'rxnConfidenceScores')
     rxnsToChange.rxnConfidenceScores=model.rxnConfidenceScores(J);
 end
+if isfield(model,'rxnDeltaG')
+    rxnsToChange.rxnDeltaG=model.rxnDeltaG(J);
+end
 if isfield(model,'pwys')
     rxnsToChange.pwys=model.pwys(J);
 end

core/checkModelStruct.m

@@ -151,6 +151,12 @@ function checkModelStruct(model,throwErrors,trimWarnings)
         dispEM(EM,throwErrors);
     end
 end
+if isfield(model,'metDeltaG')
+    if ~isnumeric(model.metDeltaG)
+        EM='The "metDeltaG" field must be a double';
+        dispEM(EM,throwErrors);
+    end
+end
 if isfield(model,'subSystems')
     for i=1:numel(model.subSystems)
         if ~iscell(model.subSystems{i,1})
@@ -189,6 +195,12 @@ function checkModelStruct(model,throwErrors,trimWarnings)
         dispEM(EM,throwErrors);
     end
 end
+if isfield(model,'rxnDeltaG')
+    if ~isnumeric(model.rxnDeltaG)
+        EM='The "rxnDeltaG" field must be a double';
+        dispEM(EM,throwErrors);
+    end
+end
 
 %Empty strings
 if isempty(model.id)
@@ -364,38 +376,38 @@ function checkModelStruct(model,throwErrors,trimWarnings)
     dispEM(EM,false,allInchis(hasMultiple),trimWarnings);
 end
 
-%Check if there are metabolites with different names but the same SMILES
-if isfield(model,'metSmiles')
-    metSmiles=containers.Map();
-    for i=1:numel(model.mets)
-        if ~isempty(model.metSmiles{i})
-            %Get existing metabolite indexes
-            if isKey(metSmiles,model.metSmiles{i})
-                existing=metSmiles(model.metSmiles{i});
-            else
-                existing=[];
-            end
-            metSmiles(model.metSmiles{i})=[existing;i];
-        end
-    end
-    
-    %Get all keys
-    allmetSmiles=keys(metSmiles);
-    
-    hasMultiple=false(numel(metSmiles),1);
-    for i=1:numel(metSmiles)
-        if numel(metSmiles(metSmiles{i}))>1
-            %Check if they all have the same name
-            if numel(unique(model.metNames(metSmiles(allmetSmiles{i}))))>1
-                hasMultiple(i)=true;
-            end
-        end
-    end
-    
-    %Print output
-    EM='The following metSmiles strings are associated to more than one unique metabolite name:';
-    dispEM(EM,false,allmetSmiles(hasMultiple),trimWarnings);
-end
+% %Check if there are metabolites with different names but the same SMILES
+% if isfield(model,'metSmiles')
+%     metSmiles=containers.Map();
+%     for i=1:numel(model.mets)
+%         if ~isempty(model.metSmiles{i})
+%             %Get existing metabolite indexes
+%             if isKey(metSmiles,model.metSmiles{i})
+%                 existing=metSmiles(model.metSmiles{i});
+%             else
+%                 existing=[];
+%             end
+%             metSmiles(model.metSmiles{i})=[existing;i];
+%         end
+%     end
+%     
+%     %Get all keys
+%     allmetSmiles=keys(metSmiles);
+%     
+%     hasMultiple=false(numel(metSmiles),1);
+%     for i=1:numel(metSmiles)
+%         if numel(metSmiles(metSmiles{i}))>1
+%             %Check if they all have the same name
+%             if numel(unique(model.metNames(metSmiles(allmetSmiles{i}))))>1
+%                 hasMultiple(i)=true;
+%             end
+%         end
+%     end
+%     
+%     %Print output
+%     EM='The following metSmiles strings are associated to more than one unique metabolite name:';
+%     dispEM(EM,false,allmetSmiles(hasMultiple),trimWarnings);
+% end
 end
 
 function I=duplicates(strings)

core/convertToIrrev.m

@@ -96,6 +96,9 @@
     if isfield(irrevModel,'rxnConfidenceScores')
         irrevModel.rxnConfidenceScores=[irrevModel.rxnConfidenceScores;irrevModel.rxnConfidenceScores(revIndexes)];
     end
+    if isfield(irrevModel,'rxnDeltaG')
+        irrevModel.rxnDeltaG=[irrevModel.rxnDeltaG;-irrevModel.rxnDeltaG(revIndexes)]; % Invert dG for reversed rxns
+    end
     if isfield(irrevModel,'rxnReferences')
         irrevModel.rxnReferences=[irrevModel.rxnReferences;irrevModel.rxnReferences(revIndexes)];
     end

core/expandModel.m

@@ -71,6 +71,9 @@
     if isfield(model,'rxnConfidenceScores')
         model.rxnConfidenceScores=[model.rxnConfidenceScores;model.rxnConfidenceScores(cpyIndices)];
     end
+    if isfield(model,'rxnDeltaG')
+        model.rxnDeltaG=[model.rxnDeltaG;model.rxnDeltaG(cpyIndices)];
+    end
 
     %now expand the more complex fields - will be filled in later
     model.rxns=[model.rxns;cell(toAdd,1)];

core/getAllowedBounds.m

@@ -1,27 +1,42 @@
-function [minFluxes, maxFluxes, exitFlags]=getAllowedBounds(model,rxns)
+function [minFluxes, maxFluxes, exitFlags]=getAllowedBounds(model,rxns,runParallel)
 % getAllowedBounds
 %   Returns the minimal and maximal fluxes through each reaction.
 %
-%   model         a model structure
-%   rxns          either a cell array of reaction IDs, a logical vector with the
-%                 same number of elements as reactions in the model, or a vector
-%                 of reaction indexes (opt, default model.rxns)
+% Input:
+%   model       a model structure
+%   rxns        either a cell array of reaction IDs, a logical vector with the
+%               same number of elements as reactions in the model, or a vector
+%               of reaction indexes (opt, default model.rxns)
+%   runParallel make use of MATLAB parallel pool to speed up calculations. Not
+%               beneficial if only a limited number of reactions are simulated.
+%               (opt, default true)
 %
-%   minFluxes     minimal allowed fluxes
-%   maxFluxes     maximal allowed fluxes
-%   exitFlags     exit flags for min/max for each of the reactions. True if
-%                 it was possible to calculate a flux
+% Output:
+%   minFluxes   minimal allowed fluxes
+%   maxFluxes   maximal allowed fluxes
+%   exitFlags   exit flags for min/max for each of the reactions. True if it was
+%               possible to calculate a flux
 %
 %   NOTE: In cases where no solution can be calculated, NaN is returned.
 %
-%   Usage: [minFluxes, maxFluxes, exitFlags]=getAllowedBounds(model,rxns)
+% Usage: [minFluxes, maxFluxes, exitFlags] = getAllowedBounds(model, rxns, runParallel)
 
 if nargin<2
     rxns=1:numel(model.rxns);
 elseif ~islogical(rxns) && ~isnumeric(rxns)
     rxns=convertCharArray(rxns);
     rxns=getIndexes(model,rxns, 'rxns');
 end
+if nargin<3
+    runParallel = true;
+end
+if runParallel
+    addonList = matlab.addons.installedAddons;
+    if ~any(strcmpi(addonList.Name,'Parallel Computing Toolbox'))
+        disp('Cannot find MATLAB Parallel Computing Toolbox, process is not parallelized.')
+        runParallel = false;
+    end
+end
 
 minFluxes=zeros(numel(rxns),1);
 maxFluxes=zeros(numel(rxns),1);
@@ -31,12 +46,9 @@
 N = numel(rxns);
 p = 1;
 h = waitbar(0, 'Please wait ...');
-
-addonList = matlab.addons.installedAddons;
-if any(strcmpi(addonList.Name,'Parallel Computing Toolbox'))
+if runParallel
     D = parallel.pool.DataQueue;
-    afterEach(D, @nUpdateWaitbar);
- 
+    afterEach(D, @nUpdateWaitbarParallel);
     parfor i=1:N
         tmpModel=model;
         tmpModel.c=c;
@@ -84,14 +96,18 @@
         else
             maxFluxes(i)=NaN;
         end
-        nUpdateWaitbar;
+        nUpdateWaitbar(p,N,h);
     end
 end
-
 close(h)
 
-    function nUpdateWaitbar(~)
-        waitbar(p/N, h);
-        p = p + 1;
-    end
+function nUpdateWaitbarParallel(~)
+waitbar(p/N, h);
+p = p + 1;
+end
+
+function nUpdateWaitbar(p,N,h)
+waitbar(p/N, h);
+p = p + 1;
+end
 end