Skip to content

Commit

Permalink
feat: prot_usage always draws from prot_pool (#379)
Browse files Browse the repository at this point in the history 
* fix: calculateFfactor can take protData input

* feat: constrainEnzConcs keep prot pool draw

addresses #375

* fix: updateProtPool obsolete

* doc: update full_ecModel for prot_usage rxns

* fix: enzymeUsage correct mention of output units

solves #376

* doc: updateGECKOdoc

* doc: update README.md with protocol change

* fix: README.md link
  • Loading branch information
@edkerk
edkerk authored May 27, 2024
1 parent 5e9ca54 commit cb54747
Show file tree
Hide file tree
Showing 10 changed files with 174 additions and 101 deletions.
3 changes: 3 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@ The **GECKO** toolbox enhances a **G**enome-scale model to account for **E**nzym

💡 In the [`GECKO/tutorials`](https://github.com/SysBioChalmers/GECKO/tree/main/tutorials) folder there are examples of how GECKO can be applied to GEMs, in either of its _full_ or _light_ forms. Each `protocol.m` contains instructions on how to reconstruct and analyze an ecModel, demonstrating how different fuctions in GECKO can be used. These two scripts complement the [protocols paper](#citation).

### Significant changes since protocol publication
- GECKO **3.2.0**: all protein usage reactions draw from the protein pool, even if they are constrained by proteomics data. This affects **Step 58** in the protocol, changing behaviour of `constrainEnzConcs` and making `updateProtPool` obsolete, `tutorials/full_ecModel/protocol.m` is updated to reflect this change. See [#357](https://github.com/SysBioChalmers/GECKO/issues/375) for more details.

_**Note:** Regarding code and model compatibility with earlier GECKO versions, see [Previous versions: GECKO 1 and 2](https://github.com/SysBioChalmers/GECKO/wiki/Previous-versions:-GECKO-1-and-2)_.

## Cite us
Expand Down
31 changes: 15 additions & 16 deletions doc/src/geckomat/limit_proteins/calculateFfactor.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,23 +121,22 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0061 uniprot = uniprotDB.ID(b(a));
0062 level(~a) = [];
0063 clear protData
0064 protData.uniprot = uniprot;
0064 protData.uniprotIDs = uniprot;
0065 protData.level = level;
0066 <span class="keyword">end</span>
0067
0068 <span class="comment">% Get MW and abundance (unit does not matter, f is fraction)</span>
0069 [~,idx] = ismember(protData.uniprot,uniprotDB.ID);
0070 protData.MW = uniprotDB.MW(idx);
0071 protData.abundance = protData.level .* protData.MW;
0072
0073 totalProt = sum(protData.abundance);
0074
0075 <span class="comment">% Get enzymes in model</span>
0076 enzymesInModel = ismember(protData.uniprot,enzymes);
0077 totalEnz = sum(protData.abundance(enzymesInModel));
0078
0079 f = totalEnz/totalProt;
0080 <span class="keyword">end</span></pre></div>
0066 <span class="comment">% Get MW and abundance (unit does not matter, f is fraction)</span>
0067 [~,idx] = ismember(protData.uniprotIDs,uniprotDB.ID);
0068 protData.MW = uniprotDB.MW(idx);
0069 protData.abundances = protData.level .* protData.MW;
0070 <span class="keyword">end</span>
0071
0072 totalProt = sum(protData.abundances);
0073
0074 <span class="comment">% Get enzymes in model</span>
0075 enzymesInModel = ismember(protData.uniprotIDs,enzymes);
0076 totalEnz = sum(protData.abundances(enzymesInModel));
0077
0078 f = totalEnz/totalProt;
0079 <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>
2 changes: 1 addition & 1 deletion doc/src/geckomat/limit_proteins/constrainEnzConcs.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0043
0044 <span class="comment">%If non-NaN in model.ec.concs, then constrain by UB</span>
0045 <span class="keyword">if</span> any(protCons)
0046 model.S(protPoolIdx, usageRxnsIdx(protCons)) = 0;
0046 <span class="comment">% model.S(protPoolIdx, usageRxnsIdx(protCons)) = 0; % Since GECKO 3.2.0</span>
0047 model.lb(usageRxnsIdx(protCons)) = -model.ec.concs(protCons);
0048 <span class="keyword">end</span>
0049 <span class="keyword">end</span></pre></div>
Expand Down
146 changes: 88 additions & 58 deletions doc/src/geckomat/limit_proteins/updateProtPool.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,13 @@ <h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> updateProtPool
Updates the protein pool to compensate for measured proteomics data (in
model.ec.concs).
Obsolete since GECKO 3.2.0, as all (measured and unmeasured) enzymes
are drawn from the protein pool. Instead, use setProtPoolSize. See
https://github.com/SysBioChalmers/GECKO/issues/375 for explanation.

Before GECKO 3.2.0: updates the protein pool to compensate for measured
proteomics data (in model.ec.concs), as only the unmeasured enzymes
draw from the protein pool.

Input:
ecModel an ecModel in GECKO 3 format (with ecModel.ec structure)
Expand All @@ -39,6 +44,7 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="
is provided, the modelAdapter value is used.
modelAdapter a loaded model adapter (Optional, will otherwise use the
default model adapter).

Output:
model an ecModel where model.ec.concs is populated with
protein concentrations
Expand All @@ -60,62 +66,86 @@ <h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^
<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function ecModel = updateProtPool(ecModel, Ptot, modelAdapter)</a>
0002 <span class="comment">% updateProtPool</span>
0003 <span class="comment">% Updates the protein pool to compensate for measured proteomics data (in</span>
0004 <span class="comment">% model.ec.concs).</span>
0005 <span class="comment">%</span>
0006 <span class="comment">% Input:</span>
0007 <span class="comment">% ecModel an ecModel in GECKO 3 format (with ecModel.ec structure)</span>
0008 <span class="comment">% Ptot total protein content in g/gDCW, overwrites the value</span>
0009 <span class="comment">% from modelAdapter. For instance, condition-specific</span>
0010 <span class="comment">% fluxData.Ptot from loadFluxData can be used. If nothing</span>
0011 <span class="comment">% is provided, the modelAdapter value is used.</span>
0012 <span class="comment">% modelAdapter a loaded model adapter (Optional, will otherwise use the</span>
0013 <span class="comment">% default model adapter).</span>
0014 <span class="comment">% Output:</span>
0015 <span class="comment">% model an ecModel where model.ec.concs is populated with</span>
0016 <span class="comment">% protein concentrations</span>
0017 <span class="comment">% Usage:</span>
0018 <span class="comment">% ecModel = updateProtPool(ecModel, Ptot, modelAdapter)</span>
0019
0020 <span class="keyword">if</span> nargin &lt; 3 || isempty(modelAdapter)
0021 modelAdapter = ModelAdapterManager.getDefault();
0022 <span class="keyword">if</span> isempty(modelAdapter)
0023 error(<span class="string">'Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.'</span>)
0024 <span class="keyword">end</span>
0025 <span class="keyword">end</span>
0026 params = modelAdapter.params;
0027
0028 <span class="keyword">if</span> nargin &lt; 2 || isempty(Ptot)
0029 Ptot = params.Ptot;
0030 disp([<span class="string">'Total protein content used: '</span> num2str(Ptot) <span class="string">' [g protein/gDw]'</span>])
0031 <span class="keyword">end</span>
0032
0033 <span class="comment">% Convert Ptot to mg/gDW if provided in g/gDCW (which is default)</span>
0034 <span class="keyword">if</span> Ptot &lt; 1
0035 Ptot = Ptot * 1000;
0036 <span class="keyword">end</span>
0037
0038 Pmeas = sum(ecModel.ec.concs,<span class="string">'omitnan'</span>);
0039 <span class="keyword">if</span> Pmeas == 0
0040 error(<span class="string">'The model has not yet been constrained with proteomics, as ecModel.ec.concs is empty.'</span>)
0041 <span class="keyword">end</span>
0042
0043 Pnew = (Ptot - Pmeas) * params.f;
0044
0045 <span class="keyword">if</span> Pnew &gt; 0
0046 PoolRxnIdx = strcmp(ecModel.rxns,<span class="string">'prot_pool_exchange'</span>);
0047 ecModel.lb(PoolRxnIdx) = -Pnew*params.sigma;
0048 sol = solveLP(ecModel);
0049 <span class="keyword">if</span> isempty(sol.x)
0050 error([<span class="string">'Estimating the remaining protein pool by subtracting the '</span> <span class="keyword">...</span>
0051 <span class="string">'sum of measured enzyme concentrations (in ecModel.ec.concs) '</span> <span class="keyword">...</span>
0052 <span class="string">'from the total protein pool (Ptot) does not yield a functional '</span> <span class="keyword">...</span>
0053 <span class="string">'model.'</span>])
0054 <span class="keyword">end</span>
0055 <span class="keyword">else</span>
0056 error(<span class="string">'The total measured protein mass exceeds the total protein content.'</span>)
0057 <span class="keyword">end</span>
0058 <span class="keyword">end</span></pre></div>
0003 <span class="comment">% Obsolete since GECKO 3.2.0, as all (measured and unmeasured) enzymes</span>
0004 <span class="comment">% are drawn from the protein pool. Instead, use setProtPoolSize. See</span>
0005 <span class="comment">% https://github.com/SysBioChalmers/GECKO/issues/375 for explanation.</span>
0006 <span class="comment">%</span>
0007 <span class="comment">% Before GECKO 3.2.0: updates the protein pool to compensate for measured</span>
0008 <span class="comment">% proteomics data (in model.ec.concs), as only the unmeasured enzymes</span>
0009 <span class="comment">% draw from the protein pool.</span>
0010 <span class="comment">%</span>
0011 <span class="comment">% Input:</span>
0012 <span class="comment">% ecModel an ecModel in GECKO 3 format (with ecModel.ec structure)</span>
0013 <span class="comment">% Ptot total protein content in g/gDCW, overwrites the value</span>
0014 <span class="comment">% from modelAdapter. For instance, condition-specific</span>
0015 <span class="comment">% fluxData.Ptot from loadFluxData can be used. If nothing</span>
0016 <span class="comment">% is provided, the modelAdapter value is used.</span>
0017 <span class="comment">% modelAdapter a loaded model adapter (Optional, will otherwise use the</span>
0018 <span class="comment">% default model adapter).</span>
0019 <span class="comment">%</span>
0020 <span class="comment">% Output:</span>
0021 <span class="comment">% model an ecModel where model.ec.concs is populated with</span>
0022 <span class="comment">% protein concentrations</span>
0023 <span class="comment">% Usage:</span>
0024 <span class="comment">% ecModel = updateProtPool(ecModel, Ptot, modelAdapter)</span>
0025
0026 <span class="comment">% Do not run from GECKO version 3.2.0 onwards. This can be recognized by</span>
0027 <span class="comment">% prot_usage reactions that are constrained by proteomics and still draw</span>
0028 <span class="comment">% from the protein pool</span>
0029 protRxns = find(startsWith(ecModel.rxns,<span class="string">'usage_prot_'</span>));
0030 poolMetIdx = find(strcmp(ecModel.mets,<span class="string">'prot_pool'</span>));
0031 <span class="comment">% Selected proteins with proteomics constraints</span>
0032 constProtRxns = ~(ecModel.lb(protRxns)==-1000);
0033 <span class="comment">% Are any still drawing from prot_pool? This is introduced in GECKO 3.2.0.</span>
0034 <span class="keyword">if</span> any(full(ecModel.S(poolMetIdx,protRxns(constProtRxns))))
0035 error([<span class="string">'In the provided ecModel, all protein usage reactions, both with '</span> <span class="keyword">...</span>
0036 <span class="string">'and without concentration constraints, draw from the protein pool. '</span> <span class="keyword">...</span>
0037 <span class="string">'This was introduced with GECKO 3.2.0. Since then, updateProtPool '</span> <span class="keyword">...</span>
0038 <span class="string">'has become obsolete, use setProtPoolSize instead to constrain the '</span> <span class="keyword">...</span>
0039 <span class="string">'total protein pool with condition-specific total protein content. See '</span><span class="keyword">...</span>
0040 <span class="string">'&lt;a href=&quot;https://github.com/SysBioChalmers/GECKO/issues/375&quot;&gt;here&lt;/a&gt; '</span> <span class="keyword">...</span>
0041 <span class="string">'for more explanation.'</span>])
0042 <span class="keyword">end</span>
0043
0044 <span class="keyword">if</span> nargin &lt; 3 || isempty(modelAdapter)
0045 modelAdapter = ModelAdapterManager.getDefault();
0046 <span class="keyword">if</span> isempty(modelAdapter)
0047 error(<span class="string">'Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.'</span>)
0048 <span class="keyword">end</span>
0049 <span class="keyword">end</span>
0050 params = modelAdapter.params;
0051
0052 <span class="keyword">if</span> nargin &lt; 2 || isempty(Ptot)
0053 Ptot = params.Ptot;
0054 disp([<span class="string">'Total protein content used: '</span> num2str(Ptot) <span class="string">' [g protein/gDw]'</span>])
0055 <span class="keyword">end</span>
0056
0057 <span class="comment">% Convert Ptot to mg/gDW if provided in g/gDCW (which is default)</span>
0058 <span class="keyword">if</span> Ptot &lt; 1
0059 Ptot = Ptot * 1000;
0060 <span class="keyword">end</span>
0061
0062 Pmeas = sum(ecModel.ec.concs,<span class="string">'omitnan'</span>);
0063 <span class="keyword">if</span> Pmeas == 0
0064 error(<span class="string">'The model has not yet been populated with proteomics, as ecModel.ec.concs is empty.'</span>)
0065 <span class="keyword">end</span>
0066
0067 Pnew = (Ptot - Pmeas) * params.f;
0068
0069 <span class="keyword">if</span> Pnew &gt; 0
0070 PoolRxnIdx = strcmp(ecModel.rxns,<span class="string">'prot_pool_exchange'</span>);
0071 ecModel.lb(PoolRxnIdx) = -Pnew*params.sigma;
0072 sol = solveLP(ecModel);
0073 <span class="keyword">if</span> isempty(sol.x)
0074 error([<span class="string">'Estimating the remaining protein pool by subtracting the '</span> <span class="keyword">...</span>
0075 <span class="string">'sum of measured enzyme concentrations (in ecModel.ec.concs) '</span> <span class="keyword">...</span>
0076 <span class="string">'from the total protein pool (Ptot) does not yield a functional '</span> <span class="keyword">...</span>
0077 <span class="string">'model.'</span>])
0078 <span class="keyword">end</span>
0079 <span class="keyword">else</span>
0080 error(<span class="string">'The total measured protein mass exceeds the total protein content.'</span>)
0081 <span class="keyword">end</span>
0082 <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>
4 changes: 2 additions & 2 deletions doc/src/geckomat/utilities/enzymeUsage.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="
<div class="fragment"><pre class="comment"> enzymeUsage
Gives enzyme usages based on a provided flux distribution, as obtained
from a full GECKO model. It can give:
1) absolute usage: the specific enzyme usage in ug/gDCW/h, which can
1) absolute usage: the specific enzyme usage in mg/gDCW, which can
be given for enzymes with- and without concentration information;
2) capacity usage: the ratio of available enzyme that is used, calcuted
by (absUsage/UB) (note that capacity usage is 0 if an enzyme
Expand Down Expand Up @@ -76,7 +76,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0002 <span class="comment">% enzymeUsage</span>
0003 <span class="comment">% Gives enzyme usages based on a provided flux distribution, as obtained</span>
0004 <span class="comment">% from a full GECKO model. It can give:</span>
0005 <span class="comment">% 1) absolute usage: the specific enzyme usage in ug/gDCW/h, which can</span>
0005 <span class="comment">% 1) absolute usage: the specific enzyme usage in mg/gDCW, which can</span>
0006 <span class="comment">% be given for enzymes with- and without concentration information;</span>
0007 <span class="comment">% 2) capacity usage: the ratio of available enzyme that is used, calcuted</span>
0008 <span class="comment">% by (absUsage/UB) (note that capacity usage is 0 if an enzyme</span>
Expand Down
17 changes: 8 additions & 9 deletions src/geckomat/limit_proteins/calculateFfactor.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,20 +61,19 @@
uniprot = uniprotDB.ID(b(a));
level(~a) = [];
clear protData
protData.uniprot = uniprot;
protData.uniprotIDs = uniprot;
protData.level = level;
% Get MW and abundance (unit does not matter, f is fraction)
[~,idx] = ismember(protData.uniprotIDs,uniprotDB.ID);
protData.MW = uniprotDB.MW(idx);
protData.abundances = protData.level .* protData.MW;
end

% Get MW and abundance (unit does not matter, f is fraction)
[~,idx] = ismember(protData.uniprot,uniprotDB.ID);
protData.MW = uniprotDB.MW(idx);
protData.abundance = protData.level .* protData.MW;

totalProt = sum(protData.abundance);
totalProt = sum(protData.abundances);

% Get enzymes in model
enzymesInModel = ismember(protData.uniprot,enzymes);
totalEnz = sum(protData.abundance(enzymesInModel));
enzymesInModel = ismember(protData.uniprotIDs,enzymes);
totalEnz = sum(protData.abundances(enzymesInModel));

f = totalEnz/totalProt;
end
2 changes: 1 addition & 1 deletion src/geckomat/limit_proteins/constrainEnzConcs.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@

%If non-NaN in model.ec.concs, then constrain by UB
if any(protCons)
model.S(protPoolIdx, usageRxnsIdx(protCons)) = 0;
% model.S(protPoolIdx, usageRxnsIdx(protCons)) = 0; % Since GECKO 3.2.0
model.lb(usageRxnsIdx(protCons)) = -model.ec.concs(protCons);
end
end
Loading

0 comments on commit cb54747

Please sign in to comment.