feat: deltaG fields (#500)

* feat: I/O YAML of deltaG fields * feat: support rxnDeltaG and metDeltaG fields * fix: rxnDeltaG and metDeltaG fields in csv files
SysBioChalmers · Jul 15, 2023 · 9d0c48a · 9d0c48a
1 parent 0e3d007
commit 9d0c48a
Show file tree

Hide file tree

Showing 40 changed files with 3,349 additions and 2,986 deletions.
diff --git a/core/addExchangeRxns.m b/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
diff --git a/core/addMets.m b/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
diff --git a/core/addRxns.m b/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

diff --git a/core/addRxnsGenesMets.m b/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

diff --git a/core/addTransport.m b/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
diff --git a/core/changeRxns.m b/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

diff --git a/core/checkModelStruct.m b/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)
@@ -363,6 +375,39 @@ function checkModelStruct(model,throwErrors,trimWarnings)
     EM='The following InChI strings are associated to more than one unique metabolite name:';
     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
 end
 
 function I=duplicates(strings)

diff --git a/core/convertToIrrev.m b/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

diff --git a/core/expandModel.m b/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)];

diff --git a/core/mergeModels.m b/core/mergeModels.m
@@ -177,6 +177,18 @@
             model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(models{i}.rxns),1)];
         end
     end
+
+    if isfield(models{i},'rxnDeltaG')
+        if isfield(model,'rxnDeltaG')
+            model.rxnDeltaG=[model.rxnDeltaG;models{i}.rxnDeltaG];
+        else
+            model.rxnDeltaG=[NaN(numel(model.rxns)-numel(models{i}.rxns),1);models{i}.rxnDeltaG];
+        end
+    else
+        if isfield(model,'rxnDeltaG')
+            model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(models{i}.rxns),1)];
+        end
+    end
 
     if isfield(models{i},'rxnComps')
         if isfield(model,'rxnComps')
@@ -359,6 +371,20 @@
                 model.metCharges=[model.metCharges;emptyMetCharge];
             end
         end
+
+        if isfield(models{i},'metDeltaG')
+            if isfield(model,'metDeltaG')
+                model.metDeltaG=[model.metDeltaG;models{i}.metDeltaG(metsToAdd)];
+            else
+                emptyMetCharge=nan(numel(model.mets)-numel(metsToAdd),1);
+                model.metDeltaG=[emptyMetCharge;models{i}.metDeltaG(metsToAdd)];
+            end
+        else
+            if isfield(model,'metDeltaG')
+                emptyMetCharge=nan(numel(metsToAdd),1);
+                model.metDeltaG=[model.metDeltaG;emptyMetCharge];
+            end
+        end
 
         if isfield(models{i},'metMiriams')
             if isfield(model,'metMiriams')

diff --git a/core/permuteModel.m b/core/permuteModel.m
@@ -76,6 +76,9 @@
         if isfield(newModel,'rxnConfidenceScores')
             newModel.rxnConfidenceScores=newModel.rxnConfidenceScores(indexes);
         end
+        if isfield(newModel,'rxnDeltaG')
+            newModel.rxnDeltaG=newModel.rxnDeltaG(indexes);
+        end
     case 'mets'
         if isfield(newModel,'mets')
             newModel.mets=newModel.mets(indexes);
@@ -113,6 +116,9 @@
         if isfield(newModel,'metCharges')
             newModel.metCharges=newModel.metCharges(indexes);
         end
+        if isfield(newModel,'metDeltaG')
+            newModel.metDeltaG=newModel.metDeltaG(indexes);
+        end
     case 'genes'
         if isfield(newModel,'genes')
             newModel.genes=newModel.genes(indexes);

diff --git a/core/predictLocalization.m b/core/predictLocalization.m
@@ -652,6 +652,9 @@
     if isfield(outModel,'rxnConfidenceScores')
         outModel.rxnConfidenceScores=[outModel.rxnConfidenceScores;outModel.rxnConfidenceScores(nER+1:nER+nRxns)];
     end
+    if isfield(outModel,'rxnDeltaG')
+        outModel.rxnDeltaG=[outModel.rxnDeltaG;outModel.rxnDeltaG(nER+1:nER+nRxns)];
+    end
     outModel.mets=[outModel.mets;strcat(outModel.mets(nEM+1:nEM+nMets),'_',GSS.compartments{i+1})];
     outModel.metNames=[outModel.metNames;outModel.metNames(nEM+1:nEM+nMets)];
     outModel.b=[outModel.b;outModel.b(nEM+1:nEM+nMets,:)];
@@ -678,6 +681,9 @@
     if isfield(outModel,'metCharges')
         outModel.metCharges=[outModel.metCharges;outModel.metCharges(nEM+1:nEM+nMets)];
     end
+    if isfield(outModel,'metDeltaG')
+        outModel.metDeltaG=[outModel.metDeltaG;outModel.metDeltaG(nEM+1:nEM+nMets)];
+    end
 end
 
 %Add the transport reactions
@@ -725,6 +731,9 @@
     if isfield(outModel,'rxnConfidenceScores')
         outModel.rxnConfidenceScores=[outModel.rxnConfidenceScores;NaN];
     end
+    if isfield(outModel,'rxnDeltaG')
+        outModel.rxnDeltaG=[outModel.rxnDeltaG;NaN];
+    end
 end
 
 %Then remove all reactions and metabolites that aren't used in the final

diff --git a/core/removeMets.m b/core/removeMets.m
@@ -92,6 +92,9 @@
     if isfield(reducedModel,'metCharges')
         reducedModel.metCharges(indexesToDelete)=[];
     end
+    if isfield(reducedModel,'metDeltaG')
+        reducedModel.metDeltaG(indexesToDelete)=[];
+    end
 end
 
 %Remove unused reactions

diff --git a/core/removeReactions.m b/core/removeReactions.m
@@ -93,6 +93,9 @@
         if isfield(reducedModel,'rxnConfidenceScores')
             reducedModel.rxnConfidenceScores(indexesToDelete,:)=[];
         end
+        if isfield(reducedModel,'rxnDeltaG')
+            reducedModel.rxnDeltaG(indexesToDelete,:)=[];
+        end
         if isfield(reducedModel,'pwys')
             reducedModel.pwys(indexesToDelete,:)=[];
         end

diff --git a/core/replaceMets.m b/core/replaceMets.m
@@ -52,6 +52,9 @@
 if isfield(model,'metCharges')
     model.metCharges(metIdx) = model.metCharges(repIdx(1));
 end
+if isfield(model,'metDeltaG')
+    model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1));
+end
 if isfield(model,'inchis')
     model.inchis(metIdx) = model.inchis(repIdx(1));
 end
@@ -98,6 +101,9 @@
     if isfield(model,'metCharges')
         model.metCharges(idxDelete) = [];
     end
+    if isfield(model,'metDeltaG')
+        model.metDeltaG(idxDelete) = [];
+    end
     if isfield(model,'inchis')
         model.inchis(idxDelete) = [];
     end

diff --git a/doc/core/addExchangeRxns.html b/doc/core/addExchangeRxns.html
@@ -153,7 +153,10 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
 0089 <span class="keyword">if</span> isfield(model,<span class="string">'rxnConfidenceScores'</span>)
 0090     model.rxnConfidenceScores=[model.rxnConfidenceScores;NaN(numel(J),1)];
 0091 <span class="keyword">end</span>
-0092 <span class="keyword">end</span></pre></div>
+0092 <span class="keyword">if</span> isfield(model,<span class="string">'rxnDeltaG'</span>)
+0093     model.rxnDeltaG=[model.rxnDeltaG;NaN(numel(J),1)];
+0094 <span class="keyword">end</span>
+0095 <span class="keyword">end</span></pre></div>
 <hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
 </body>
 </html>