Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix/unbalanced_rxns #222

Merged
merged 16 commits into from
Jun 12, 2020
Merged
Show file tree
Hide file tree
Changes from 14 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
843 changes: 843 additions & 0 deletions ComplementaryData/modelCuration/modMetsandSmatrix.tsv

Large diffs are not rendered by default.

68 changes: 68 additions & 0 deletions ComplementaryScripts/modelCuration/changerxn.m
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
102 changes: 102 additions & 0 deletions ComplementaryScripts/modelCuration/checkMassChargeBalance.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
% 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 issue causing imbalance(s),
% individual charges of each metabolite in the reaction, sum of all
% charges in the reaction and elemental difference will be
% displayed in MassChargeresults
%
% NOTE: getElementalBalance.m is a function from RAVEN
%
% modified from Feiran Li's script 'checkBalanceforSce.m'
%
% Cheng Wei Quan (Eiden), 2020-05-06
eiden309 marked this conversation as resolved.
Show resolved Hide resolved

function MassChargeresults = checkMassChargeBalance(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