Fix/unbalanced_rxns

ComplementaryScripts/modelCuration/changerxn.m

@@ -0,0 +1,68 @@
+function model = changerxn(model,rxnID,rxnformula,grRule)
+% changerxn
+%   Change the metabolites in the reaction
+%
+%   model           COBRA model structure
+%   rxnID           ID of the reaction i.e. identifier in model.rxns
+%   rxnFormula      Reaction formula in string format 
+%                   (A [comp] + B [comp] -> C [comp])
+%                   For metabolites with space in model.metNames e.g. 
+%                   1-pyrroline-3-hydroxy-5-carboxylic acid, use symbol '&'
+%                   to substitute space ' ' i.e. 1-pyrroline-3-hydroxy-5-carboxylic&acid
+%   grRule          Gene-reaction rule in boolean format (and/or allowed)
+%                   (opt, default model.grRule from RAVEN model structure)               
+%
+%   Note: model.mets is generated for new metabolites with prefix 's_'
+%
+%   Usage: model = changerxn(model,rxnID,rxnformula,grRule)
+%
+% Feiran Li, 2020-05-24
+% Cheng Wei Quan (Eiden), added documentation, 2020-05-24
+
+[~,idx] = ismember(rxnID,model.rxns); 
+if nargin < 4
+    if isfield(model,'grRules')
+        grRule = model.grRules{idx};
+    else
+        model1 = ravenCobraWrapper(model);
+         grRule = model1.grRules{idx};       
+    end
+end
+
+rxnformula = strrep(rxnformula,' [','[');
+[metaboliteList, stoichCoeffList, revFlag] = parseRxnFormula(rxnformula);
+metaboliteList = strrep(metaboliteList,'[',' [');
+metaboliteList = strrep(metaboliteList,'&',' ');
+comps = split(metaboliteList', ' [');
+comps = comps(:,2);
+comps = strrep(comps,']','');
+CONValldata = cat(2,model.compNames,model.comps);
+[~,b] = ismember(comps,CONValldata(:,1));
+comps = CONValldata(b,2);
+
+%mapping mets to model.metnames, get s_ index for new mets
+for j = 1:length(metaboliteList)
+    [~,metindex] = ismember(metaboliteList(j),model.metNames);
+    if metindex ~= 0
+        mets(j) = model.mets(metindex);
+    elseif metindex == 0
+        cd ../otherChanges/
+        newID = getNewIndex(model.mets);
+        cd ../modelCuration/
+        mets(j) = strcat('s_',newID,'[',comps(j),']');
+        model = addMetabolite(model,char(mets(j)), ...
+                            'metName',metaboliteList{j});
+    end
+end
+
+
+[model, rxnIDexists] = addReaction(model,...
+                                    rxnID,...
+                                    'reactionName', model.rxnNames{idx},...
+                                    'metaboliteList',mets,...
+                                    'stoichCoeffList',stoichCoeffList,...
+                                    'reversible',revFlag,...
+                                    'geneRule',grRule,...
+                                    'checkDuplicate',1);
+
+end

ComplementaryScripts/modelCuration/checkRxnBalance.m

@@ -0,0 +1,104 @@
+% This function is to identify charge and mass balance for reactions
+%
+% Input: model and reaction identifier(s) in model.rxns (e.g. 'r_0001')
+% Output: MassChargeresults
+%         For exchange reactions:
+%         reaction identifier(s) and 'exchange' will be displayed in
+%         MassChargeresults
+%
+%         For balanced reactions:
+%         reaction identifier(s) and 'pass' will be displayed in
+%         MassChargeresults
+%
+%         For reaction identifier(s) which does not match model.rxns: 
+%         'rxnID not found' and 'manual check required' will be displayed
+%         in MassChargeresults
+%
+%         For reactions with elemental and/or charge imbalance(s):  
+%         reaction identifier(s), specific imbalance issue(s) to check, 
+%         individual charges of each metabolite in the reaction (if there 
+%         is charge difference), charge difference and elemental difference
+%         will be displayed in MassChargeresults
+%
+% NOTE: getElementalBalance.m is a function from RAVEN
+%
+% Usage: MassChargeresults = checkRxnBalance(model,rxn)
+%
+% modified from Feiran Li's script 'checkBalanceforSce.m'
+%
+% Cheng Wei Quan (Eiden), 2020-05-06
+
+function MassChargeresults = checkRxnBalance(model,rxn)
+exchangeRxns = findExcRxns(model);
+MassChargeresults{length(rxn),5} = [];
+for i = 1:length(rxn)
+    [~,rxnID] = ismember(rxn(i),model.rxns);
+    if rxnID ~= 0
+        if exchangeRxns(rxnID) == 1
+            arrayidx = find(cellfun('isempty', MassChargeresults),1);
+            MassChargeresults(arrayidx,1) = model.rxns(rxnID);
+            MassChargeresults(arrayidx,2) = {'exchange'};
+        else
+            %check mass balance
+            balanceStructure=getElementalBalance(model,rxn(i));
+            unbalancedidx = find(balanceStructure.balanceStatus==0);
+            missingidx = find(balanceStructure.balanceStatus==-1);
+            idx = union(unbalancedidx,missingidx);
+            for j=1:length(idx)
+                dif = balanceStructure.leftComp(idx(j),:) - balanceStructure.rightComp(idx(j),:);
+                mets = find(~dif==0);
+                difference = '';
+                for k=1:length(mets)
+                    difference = strcat(difference,balanceStructure.elements.abbrevs(mets(k)));
+                    difference = strcat(difference,num2str(dif(mets(k))));
+                end
+                out = difference;
+            end
+            %check charge balance
+            mets=find(model.S(:,rxnID));
+            coef = model.metCharges(mets);
+            if length(mets) == length(coef)
+                indvCharge = model.S(mets,rxnID).*coef;
+                balanceCharge = sum(model.S(mets,rxnID).*coef);
+            else
+                balanceCharge = -1;
+            end
+            if balanceStructure.balanceStatus == 1 && balanceCharge == 0
+                arrayidx = find(cellfun('isempty', MassChargeresults),1);
+                MassChargeresults(arrayidx,1) = model.rxns(rxnID);
+                MassChargeresults(arrayidx,2) = {'pass'};
+            elseif balanceStructure.balanceStatus == 1 && balanceCharge ~= 0
+                arrayidx = find(cellfun('isempty', MassChargeresults),1);
+                MassChargeresults(arrayidx,1) = model.rxns(rxnID);
+                MassChargeresults(arrayidx,2) = {'checkMetCharge'};
+                MassChargeresults(arrayidx,3) = {indvCharge};
+                MassChargeresults(arrayidx,4) = {balanceCharge};
+            elseif balanceStructure.balanceStatus == 0 && balanceCharge == 0
+                arrayidx = find(cellfun('isempty', MassChargeresults),1);
+                MassChargeresults(arrayidx,1) = model.rxns(rxnID);
+                MassChargeresults(arrayidx,2) = {'checkMetFormula'};
+                MassChargeresults(arrayidx,5) = out;
+            elseif balanceStructure.balanceStatus == 0 && balanceCharge ~= 0
+                arrayidx = find(cellfun('isempty', MassChargeresults),1);
+                MassChargeresults(arrayidx,1) = model.rxns(rxnID);
+                MassChargeresults(arrayidx,2) = {'checkMetChargeAndFormula'};
+                MassChargeresults(arrayidx,3) = {indvCharge};
+                MassChargeresults(arrayidx,4) = {balanceCharge};
+                MassChargeresults(arrayidx,5) = out;
+            elseif balanceStructure.balanceStatus == -1
+                arrayidx = find(cellfun('isempty', MassChargeresults),1);
+                MassChargeresults(arrayidx,1) = model.rxns(rxnID);
+                MassChargeresults(arrayidx,2) = {'unbalanced due to missing information'};
+                MassChargeresults(arrayidx,3) = {indvCharge};
+                MassChargeresults(arrayidx,4) = {balanceCharge};
+                MassChargeresults(arrayidx,5) = {'-'};
+            end
+        end
+    else
+        arrayidx = find(cellfun('isempty', MassChargeresults),1);
+        MassChargeresults(arrayidx,1) = {'rxnID not found'};
+        MassChargeresults(arrayidx,2) = {'manual check required'};
+    end
+end
+
+end