Skip to content

Commit

Permalink
Merge pull request #222 from SysBioChalmers/fix/unbalanced_rxns
Browse files Browse the repository at this point in the history 
Fix/unbalanced_rxns
  • Loading branch information
@BenjaSanchez
BenjaSanchez authored Jun 12, 2020
2 parents 984dca1 + 2ef9281 commit 39ada1b
Show file tree
Hide file tree
Showing 12 changed files with 6,295 additions and 2,267 deletions.
843 changes: 843 additions & 0 deletions ComplementaryData/modelCuration/modMetsandSmatrix.tsv
View file

Large diffs are not rendered by default.

68 changes: 68 additions & 0 deletions ComplementaryScripts/modelCuration/changerxn.m
View file
Original file line number Diff line number Diff line change
@@ -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
104 changes: 104 additions & 0 deletions ComplementaryScripts/modelCuration/checkRxnBalance.m
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 39ada1b

Please sign in to comment.