Skip to content

Commit

Permalink
fix: getSubsetEcModel check for reaction presence (#359)
Browse files Browse the repository at this point in the history
  • Loading branch information
@edkerk
edkerk authored Dec 2, 2023
1 parent 9546b87 commit 074078c
Show file tree
Hide file tree
Showing 2 changed files with 79 additions and 45 deletions.
110 changes: 65 additions & 45 deletions doc/src/geckomat/utilities/getSubsetEcModel.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,12 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="
Generate a context-specific ecModel (strain/cell-line/tissue) by
mapping reactions and genes from a conventional context-specific GEM to
a general ecModel (bigEcModel), to yield a context-specific ecModel.

Both models (bigEcModel and the smallGEM) should have derived from the
same starting model. If an ecModel with human-GEM 1.15.0 is made
(yielding the bigEcModel), then human-GEM 1.15.0 should also be the
basis from which the context-specific GEM (smallGEM) was generated by
for instance tINIT.

Input:
bigEcModel enzyme-constrained version of conventional model of
Expand Down Expand Up @@ -62,53 +68,67 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0004 <span class="comment">% mapping reactions and genes from a conventional context-specific GEM to</span>
0005 <span class="comment">% a general ecModel (bigEcModel), to yield a context-specific ecModel.</span>
0006 <span class="comment">%</span>
0007 <span class="comment">% Input:</span>
0008 <span class="comment">% bigEcModel enzyme-constrained version of conventional model of</span>
0009 <span class="comment">% metabolism for a given species</span>
0010 <span class="comment">% smallGEM Reduced model (subset of the general model) for a</span>
0011 <span class="comment">% specific strain (microbes) or cell-line/tissues (mammals)</span>
0007 <span class="comment">% Both models (bigEcModel and the smallGEM) should have derived from the</span>
0008 <span class="comment">% same starting model. If an ecModel with human-GEM 1.15.0 is made</span>
0009 <span class="comment">% (yielding the bigEcModel), then human-GEM 1.15.0 should also be the</span>
0010 <span class="comment">% basis from which the context-specific GEM (smallGEM) was generated by</span>
0011 <span class="comment">% for instance tINIT.</span>
0012 <span class="comment">%</span>
0013 <span class="comment">% Output:</span>
0014 <span class="comment">% smallEcModel Enzyme-constrained version of the context-specific model</span>
0015 <span class="comment">%</span>
0016 <span class="comment">% Usage: smallEcModel = getSubsetEcModel(smallGEM,bigEcModel)</span>
0017
0018 <span class="comment">% Check if original bigEcModel contains context-dependent protein constraints</span>
0019 <span class="keyword">if</span> any(bigEcModel.lb(startsWith(bigEcModel.rxns,<span class="string">'usage_prot_'</span>)) ~= -1000)
0020 printOrange([<span class="string">'WARNING: The bigEcModel is constraint by protein concentrations that are\n'</span> <span class="keyword">...</span>
0021 <span class="string">'likely not relevant in the constructed smallEcModel.\n'</span>])
0022 <span class="keyword">end</span>
0013 <span class="comment">% Input:</span>
0014 <span class="comment">% bigEcModel enzyme-constrained version of conventional model of</span>
0015 <span class="comment">% metabolism for a given species</span>
0016 <span class="comment">% smallGEM Reduced model (subset of the general model) for a</span>
0017 <span class="comment">% specific strain (microbes) or cell-line/tissues (mammals)</span>
0018 <span class="comment">%</span>
0019 <span class="comment">% Output:</span>
0020 <span class="comment">% smallEcModel Enzyme-constrained version of the context-specific model</span>
0021 <span class="comment">%</span>
0022 <span class="comment">% Usage: smallEcModel = getSubsetEcModel(smallGEM,bigEcModel)</span>
0023
0024 <span class="comment">% Remove genes (and associated reactions) that are absent in smallGEM</span>
0025 genesToRemove = setdiff(bigEcModel.genes,smallGEM.genes);
0026 smallEcModel = removeGenes(bigEcModel,genesToRemove,true,true,false);
0027
0028 <span class="comment">% Remove genes from ec-structure</span>
0029 enzToRemove = ismember(smallEcModel.ec.genes,genesToRemove);
0030 smallEcModel.ec.genes(enzToRemove) = [];
0031 smallEcModel.ec.enzymes(enzToRemove) = [];
0032 smallEcModel.ec.concs(enzToRemove) = [];
0033 smallEcModel.ec.mw(enzToRemove) = [];
0034 smallEcModel.ec.sequence(enzToRemove) = [];
0035 smallEcModel.ec.rxnEnzMat(:,enzToRemove) = [];
0036
0037 <span class="comment">% Remove any reaction (except prot_ associated) that is absent from smallGEM</span>
0038 trimRxns = regexprep(smallEcModel.rxns,<span class="string">'_REV|_EXP_\d+'</span>,<span class="string">''</span>);
0039 protRxns = startsWith(smallEcModel.rxns,{<span class="string">'usage_prot_'</span>,<span class="string">'prot_pool_exchange'</span>});
0040 keepRxns = ismember(trimRxns,smallGEM.rxns) | protRxns;
0041 smallEcModel = removeReactions(smallEcModel,~keepRxns, true, true, true);
0042
0043 <span class="comment">% Remove reactions from ec-structure</span>
0044 ecRxnsToRemove = ~ismember(smallEcModel.ec.rxns,smallEcModel.rxns);
0045
0046 smallEcModel.ec.rxns(ecRxnsToRemove) = [];
0047 smallEcModel.ec.kcat(ecRxnsToRemove) = [];
0048 smallEcModel.ec.source(ecRxnsToRemove) = [];
0049 smallEcModel.ec.notes(ecRxnsToRemove) = [];
0050 smallEcModel.ec.eccodes(ecRxnsToRemove) = [];
0051 smallEcModel.ec.rxnEnzMat(ecRxnsToRemove,:) = [];
0052 <span class="keyword">end</span>
0053</pre></div>
0024 <span class="comment">% Check if models were derived from the same starting GEM</span>
0025 rxnsDiff = find(~ismember(smallGEM.rxns,bigEcModel.rxns));
0026 <span class="keyword">if</span> numel(rxnsDiff) &gt; 0
0027 dispEM([<span class="string">'While both models should have derived from the same starting '</span>,<span class="keyword">...</span>
0028 <span class="string">'GEM, the following reactions from smallGEM could not be found '</span>,<span class="keyword">...</span>
0029 <span class="string">'in bigEcModel:'</span>],true,smallGEM.rxns(rxnsDiff),false);
0030 <span class="keyword">end</span>
0031
0032 <span class="comment">% Check if original bigEcModel contains context-dependent protein constraints</span>
0033 <span class="keyword">if</span> any(bigEcModel.lb(startsWith(bigEcModel.rxns,<span class="string">'usage_prot_'</span>)) ~= -1000)
0034 printOrange([<span class="string">'WARNING: The bigEcModel is constraint by protein concentrations that are\n'</span> <span class="keyword">...</span>
0035 <span class="string">'likely not relevant in the constructed smallEcModel.\n'</span>])
0036 <span class="keyword">end</span>
0037
0038 <span class="comment">% Remove genes (and associated reactions) that are absent in smallGEM</span>
0039 genesToRemove = setdiff(bigEcModel.genes,smallGEM.genes);
0040 smallEcModel = removeGenes(bigEcModel,genesToRemove,true,true,false);
0041
0042 <span class="comment">% Remove genes from ec-structure</span>
0043 enzToRemove = ismember(smallEcModel.ec.genes,genesToRemove);
0044 smallEcModel.ec.genes(enzToRemove) = [];
0045 smallEcModel.ec.enzymes(enzToRemove) = [];
0046 smallEcModel.ec.concs(enzToRemove) = [];
0047 smallEcModel.ec.mw(enzToRemove) = [];
0048 smallEcModel.ec.sequence(enzToRemove) = [];
0049 smallEcModel.ec.rxnEnzMat(:,enzToRemove) = [];
0050
0051 <span class="comment">% Remove any reaction (except prot_ associated) that is absent from smallGEM</span>
0052 trimRxns = regexprep(smallEcModel.rxns,<span class="string">'_REV|_EXP_\d+'</span>,<span class="string">''</span>);
0053 protRxns = startsWith(smallEcModel.rxns,{<span class="string">'usage_prot_'</span>,<span class="string">'prot_pool_exchange'</span>});
0054 keepRxns = ismember(trimRxns,smallGEM.rxns) | protRxns;
0055 smallEcModel = removeReactions(smallEcModel,~keepRxns, true, true, true);
0056
0057 <span class="comment">% Remove reactions from ec-structure</span>
0058 ecRxnsToRemove = ~ismember(smallEcModel.ec.rxns,smallEcModel.rxns);
0059
0060 smallEcModel.ec.rxns(ecRxnsToRemove) = [];
0061 smallEcModel.ec.kcat(ecRxnsToRemove) = [];
0062 smallEcModel.ec.source(ecRxnsToRemove) = [];
0063 smallEcModel.ec.notes(ecRxnsToRemove) = [];
0064 smallEcModel.ec.eccodes(ecRxnsToRemove) = [];
0065 smallEcModel.ec.rxnEnzMat(ecRxnsToRemove,:) = [];
0066 <span class="keyword">end</span>
0067</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>
14 changes: 14 additions & 0 deletions src/geckomat/utilities/getSubsetEcModel.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,12 @@
% Generate a context-specific ecModel (strain/cell-line/tissue) by
% mapping reactions and genes from a conventional context-specific GEM to
% a general ecModel (bigEcModel), to yield a context-specific ecModel.
%
% Both models (bigEcModel and the smallGEM) should have derived from the
% same starting model. If an ecModel with human-GEM 1.15.0 is made
% (yielding the bigEcModel), then human-GEM 1.15.0 should also be the
% basis from which the context-specific GEM (smallGEM) was generated by
% for instance tINIT.
%
% Input:
% bigEcModel enzyme-constrained version of conventional model of
Expand All @@ -15,6 +21,14 @@
%
% Usage: smallEcModel = getSubsetEcModel(smallGEM,bigEcModel)

% Check if models were derived from the same starting GEM
rxnsDiff = find(~ismember(smallGEM.rxns,bigEcModel.rxns));
if numel(rxnsDiff) > 0
dispEM(['While both models should have derived from the same starting ',...
'GEM, the following reactions from smallGEM could not be found ',...
'in bigEcModel:'],true,smallGEM.rxns(rxnsDiff),false);
end

% Check if original bigEcModel contains context-dependent protein constraints
if any(bigEcModel.lb(startsWith(bigEcModel.rxns,'usage_prot_')) ~= -1000)
printOrange(['WARNING: The bigEcModel is constraint by protein concentrations that are\n' ...
Expand Down

0 comments on commit 074078c

Please sign in to comment.