diff --git a/GECKOInstaller.m b/GECKOInstaller.m index 434b96952..2811abc84 100644 --- a/GECKOInstaller.m +++ b/GECKOInstaller.m @@ -28,6 +28,7 @@ end addpath(paths); savepath; + GECKOInstaller.checkGECKOversion; end function uninstall sourceDir = fileparts(which(mfilename)); @@ -59,6 +60,7 @@ pathsLeft = splitPaths(1,okPaths); newPaths = char(join(pathsLeft, pathSep)); end + function checkRAVENversion(minmVer) try currVer = checkInstallation('versionOnly'); @@ -81,6 +83,39 @@ function checkRAVENversion(minmVer) end end + + function checkGECKOversion + sourceDir = fileparts(which(mfilename)); + hasGit=isfolder(fullfile(sourceDir,'.git')); + + if exist(fullfile(sourceDir,'version.txt'), 'file') == 2 + currVer = fgetl(fopen(fullfile(sourceDir,'version.txt'))); + fclose('all'); + fprintf('GECKO version %s installed',currVer) + try + newVer=strtrim(webread('https://raw.githubusercontent.com/SysBioChalmers/GECKO/main/version.txt')); + newVerNum=str2double(strsplit(newVer,'.')); + currVerNum=str2double(strsplit(currVer,'.')); + for i=1:3 + if currVerNum(i)here\n'); + end + break + elseif i==3 + fprintf('\n'); + end + end + catch + fprintf('\n'); + end + else + fprintf('GECKO installed, unknown version (cannot find version.txt)\n') + end + end end properties (Constant) diff --git a/README.md b/README.md index c0fee9c1f..18cf48537 100644 --- a/README.md +++ b/README.md @@ -1,72 +1,26 @@ ![Current version](https://badge.fury.io/gh/sysbiochalmers%2Fgecko.svg) -[![Gitter](https://badges.gitter.im/SysBioChalmers/GECKO.svg)](https://gitter.im/SysBioChalmers/GECKO) +[![GitHub Discussions](https://img.shields.io/github/discussions-search?query=repo%3Asysbiochalmers%2Fgecko&label=GitHub%20Discussions)](https://github.com/SysBioChalmers/GECKO/discussions) [![Zenodo](https://zenodo.org/badge/DOI/10.5281/zenodo.7699818.svg)](https://doi.org/10.5281/zenodo.7699818) ## About GECKO 3 The **GECKO** toolbox enhances a **G**enome-scale model to account for **E**nzyme **C**onstraints, using **K**inetics and **O**mics. The resulting enzyme-constrained model (**ecModel**) can be used to perform simulations where enzyme allocation is either drawn from a total protein pool, or constrained by measured protein levels from proteomics data. -💡 In the GECKO folder, `protocol.m` contains instructions on how to reconstruct and analyze an ecModel for _S. cerevisiae_. This demonstrates how many of GECKO's functions can be used. -💡 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. The source code documentation is also available -[online](http://sysbiochalmers.github.io/GECKO/doc/). +💡 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). -_**Note:** Regarding code and model compatibility with earlier GECKO versions, see [Previous versions: GECKO 1 and 2](#previous-versions-gecko-1-and-2)_. +_**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)_. -### Citation -- A GECKO 3 publication is currently under consideration, citation information will appear here in due course. -- For GECKO release 2, please cite [Domenzain et al. (2022) doi:10.1038/s41467-022-31421-1](https://doi.org/10.1038/s41467-022-31421-1). -- For GECKO release 1, please cite [Sánchez et al. (2017) doi:10.15252/msb.20167411](https://doi.org/10.15252/msb.20167411). +## Documentation +**Installation instructions** and other information is available in the **[Wiki](https://github.com/SysBioChalmers/GECKO/wiki)**. The source code documentation is also available +[online](http://sysbiochalmers.github.io/GECKO/doc/). Use [GitHub Discussions](https://github.com/SysBioChalmers/GECKO/discussions) for support, to ask questions or leave comments. -### Getting started +## Citation -#### Required software - -- MATLAB version 2019b or later, no additional MathWorks toolboxes are required. -- [RAVEN Toolbox](https://github.com/SysBioChalmers/RAVEN). The RAVEN Toolbox Wiki contains [installation instructions for both RAVEN](https://github.com/SysBioChalmers/RAVEN/wiki/Installation) and [Gurobi](https://github.com/SysBioChalmers/RAVEN/wiki/Installation#solvers). Briefly, RAVEN is either downloaded via `git clone`, as ZIP-archive from GitHub, or installed as a [MATLAB AddOn](https://se.mathworks.com/matlabcentral/fileexchange/112330-raven-toolbox). After finishing all installation instructions, the user should run installation checks in MATLAB with: `checkInstallation`. -- [Gurobi Optimizer](https://www.gurobi.com/solutions/gurobi-optimizer/) is recommended for simulations (free academic license available). Alternatively, the open-source [GNU Linear Programming Kit](https://www.gnu.org/software/glpk/) (distributed with RAVEN) or SoPlex as part of the [SCIP Optimization Suite](https://scipopt.org/) can be used. -- [Docker](https://www.docker.com/) for running DLKcat. Installation instructions are available at https://docs.docker.com/get-docker . - -#### Installation - -- The preferred way to download GECKO is via git clone: - -``` -git clone --depth=1 https://github.com/SysBioChalmers/GECKO -``` - -- Alternatively, [a ZIP-archive can be directly downloaded from GitHub](https://github.com/SysBioChalmers/GECKO/releases). The ZIP-archive should be extracted to a disk location where the user has read- and write-access rights. - -- After `git clone` or extracting the ZIP-archive, the user should navigate in MATLAB to the GECKO folder. GECKO can then be installed with the command that adds GECKO (sub-)folders to the MATLAB path:: - -``` - cd('C:\path\to\GECKO') % Modify to match GECKO folder and operating system - GECKOInstaller.install -``` - -- If desired, a removal command is available as:: - -``` - GECKOInstaller.uninstall -``` - -All set! 🚀 - -## Previous versions: GECKO 1 and 2 - -Due to significant refactoring of the code, GECKO version 3 is largely not backwards compatible with earlier GECKO versions. - -- Most notably, GECKO 3 ecModels have an `.ec` structure containing all enzyme and kcat information. -- In addition, in GECKO 3 enzymes are incorporated in the S-matrix as MW/kcat, while in GECKO 1 and 2 this was 1/kcat (where the MW was instead considered in the protein exchange reactions). -- GECKO 3 ecModels can be stored in YAML file format that retains all model content. -- Most functions in GECKO 3 do not work on ecModels generated with GECKO versions 1 or 2. -- ecModels generated in GECKO 3 do not work with functions from GECKO versions 1 or 2. -- At this moment, there are no Python functions to work with GECKO 3 formatted ecModels. -- The last GECKO 2 release (2.0.3) is available [here](https://github.com/SysBioChalmers/GECKO/releases/tag/v2.0.3). -- The `gecko2` branch remains available for any potential fixes. +- The GECKO 3 publication is in press at Nature Protocols. Once published, citation information will appear here. ## Contributing diff --git a/doc/index.html b/doc/index.html index b73088c93..8a709fdb6 100644 --- a/doc/index.html +++ b/doc/index.html @@ -19,23 +19,22 @@

Matlab Directories

Matlab Files found in these Directories

- - - - - - - - - - - - - - - - -
ModelAdapter ecFSEOF getReactionsFromEnzyme runDLKcat
ModelAdapterManager ecFVA getStandardKcat saveEcModel
abc_max enzymeUsage getSubsetEcModel selectKcatValue
adapterTemplate fillEnzConcs getcarbonnum sensitivityTuning
addCarbonNum findECInDB getrSample setKcatForReactions
addNewRxnsToEC findGECKOroot loadBRENDAdata setProtPoolSize
anaerobicModel findMaxValue loadConventionalGEM sigmaFitter
applyComplexData findMetSmiles loadDatabases startGECKOproject
applyCustomKcats flexibilizeEnzConcs loadEcModel truncateValues
applyKcatConstraints fuzzyKcatMatching loadFluxData updateGECKOdoc
bayesianSensitivityTuning getComplexData loadProtData updateProtPool
calculateFfactor getConcControlCoeffs makeEcModel updateprior
calculateMW getECfromDatabase mapRxnsToConv writeDLKcatInput
changeMedia getECfromGEM mergeDLKcatAndFuzzyKcats
constrainEnzConcs getECstring plotEcFVA
constrainFluxData getFluxTarget readDLKcatOutput
copyECtoGEM getKcatAcrossIsozymes reportEnzymeUsage
+ ModelAdapter copyECtoGEM getKcatAcrossIsozymes readDLKcatOutput + ModelAdapterManager ecFSEOF getReactionsFromEnzyme reportEnzymeUsage + abc_max ecFVA getStandardKcat runDLKcat + adapterTemplate enzymeUsage getSubsetEcModel saveEcModel + addCarbonNum fillEnzConcs getcarbonnum selectKcatValue + addNewRxnsToEC findECInDB getrSample sensitivityTuning + anaerobicModel findGECKOroot loadBRENDAdata setKcatForReactions + applyComplexData findMaxValue loadConventionalGEM setProtPoolSize + applyCustomKcats findMetSmiles loadDatabases sigmaFitter + applyKcatConstraints flexibilizeEnzConcs loadEcModel startGECKOproject + bayesianSensitivityTuning fuzzyKcatMatching loadFluxData truncateValues + calculateFfactor getComplexData loadProtData updateGECKOdoc + calculateMW getConcControlCoeffs makeEcModel updateProtPool + changeMedia getECfromDatabase mapRxnsToConv updateprior + constrainEnzConcs getECfromGEM mergeDLKcatAndFuzzyKcats writeDLKcatInput + constrainFluxData getECstring plotEcFVA
Generated by m2html © 2005
diff --git a/doc/src/geckomat/change_model/makeEcModel.html b/doc/src/geckomat/change_model/makeEcModel.html index cf51d3bcb..18aeab8e1 100644 --- a/doc/src/geckomat/change_model/makeEcModel.html +++ b/doc/src/geckomat/change_model/makeEcModel.html @@ -463,7 +463,7 @@

SOURCE CODE ^%9: Add proteins as pseudometabolites 0328 if ~geckoLight -0329 [proteinMets.mets, uniprotSortId] = sort(ec.enzymes); +0329 [proteinMets.mets, uniprotSortId] = unique(ec.enzymes); 0330 proteinMets.mets = strcat('prot_',proteinMets.mets); 0331 proteinMets.metNames = proteinMets.mets; 0332 proteinMets.compartments = compartmentID; diff --git a/doc/src/geckomat/utilities/ecFSEOF.html b/doc/src/geckomat/utilities/ecFSEOF.html index eb327cf71..c6aaab129 100644 --- a/doc/src/geckomat/utilities/ecFSEOF.html +++ b/doc/src/geckomat/utilities/ecFSEOF.html @@ -24,7 +24,7 @@

PURPOSE ^ecFSEOF

SYNOPSIS ^

-
function FC = ecFSEOF(ecModel,targetRxn,csRxn,alphaLims,nSteps,filePath,filterG,modelAdapter)
+
function fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,modelAdapter)

DESCRIPTION ^

 ecFSEOF
@@ -32,52 +32,52 @@ 

DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

This function calls: + This function is called by:
@@ -86,54 +86,54 @@

CROSS-REFERENCE INFORMATION ^
 
 
 <h2><a name=SOURCE CODE ^

-
0001 function FC = ecFSEOF(ecModel,targetRxn,csRxn,alphaLims,nSteps,filePath,filterG,modelAdapter)
+
0001 function fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,modelAdapter)
 0002 % ecFSEOF
 0003 %   Function that runs Flux-Scanning with Enforced Objective Function (FSEOF)
 0004 %   for a specified production target.
 0005 %
 0006 % Input:
-0007 %   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
-0008 %   targetRxn       rxn ID for the production target reaction, a exchange
-0009 %                   reaction is recommended.
-0010 %   csRxn           rxn ID for the main carbon source uptake reaction.
-0011 %   alphaLims       vector of Minimum and maximum biomass scalling factors for
-0012 %                   enforced objective limits (e.g. [0.5 0.75]). Max value: 1.
-0013 %                   Max value recomended 0.9 (Optional, default [0.5 0.75])
-0014 %   nSteps          number of steps for suboptimal objective in FSEOF.
-0015 %                   (Optional, default 16)
-0016 %   filePath        file path for results output. It will store two files:
-0017 %                   - at the genes level, ecFSEOF_genes.tsv
-0018 %                   - at the reactions level, ecFSEOF_rxns.tsv
-0019 %                   (Optional, default in the 'output' sub-folder taken from
-0020 %                   modelAdapter, e.g. output/ecFSEOF_rxns.tsv)
-0021 %   filterG         logical value. TRUE if genes K_scores results should be
-0022 %                   filtered according to the alpha vector distribution
-0023 %                   (Optional, default true)
-0024 %   modelAdapter    a loaded model adapter. (Optional, will otherwise use
-0025 %                   the default model adapter)
-0026 %
-0027 % Output:
-0028 %   FC           structure with all results. Contains the following fields:
-0029 %                - flux_WT:   flux distribution for the WT strain (100% of
-0030 %                             carbon towards growth).
-0031 %                - alpha:     scalling factors used for enforced objetive
-0032 %                             limits
-0033 %                - v_matrix:  fluxes for each reaction in rxnsTable.rxns
-0034 %                             and each alpha.
-0035 %                - k_matrix:  fold changes for each reaction in
-0036 %                             rxnsTable.rxns and each alpha.
-0037 %                - rxnsTable: a list with all reactions with fluxes that
-0038 %                             change consistently as target production
-0039 %                             increases.
-0040 %                             Contains: ID, name, k_score, gene rule, equation
-0041 %                - geneTable: a list with all selected targets that
-0042 %                             increase production.
-0043 %                             Contains: gene, shortName, k_score
+0007 %   model          an ecModel in GECKO 3 format (with ecModel.ec structure).
+0008 %   prodTargetRxn  rxn ID for the production target reaction, a exchange
+0009 %                  reaction is recommended.
+0010 %   csRxn          rxn ID for the main carbon source uptake reaction.
+0011 %   nSteps         number of steps for suboptimal objective in FSEOF.
+0012 %                  (Optional, default 16)
+0013 %   outputFile     bolean option to save results in a file. (Optional,
+0014 %                  default false)
+0015 %   filePath       file path for results output. It will store two files:
+0016 %                  - at the genes level, ecFSEOF_genes.tsv
+0017 %                  - at the reactions level, ecFSEOF_rxns.tsv
+0018 %                  (Optional, default in the 'output' sub-folder taken from
+0019 %                  modelAdapter, e.g. output/ecFSEOF_rxns.tsv)
+0020 %   modelAdapter   a loaded model adapter. (Optional, will otherwise use
+0021 %                  the default model adapter)
+0022 %
+0023 % Output:
+0024 %   fseof   an structure with all results. Contains the following fields:
+0025 %           - alpha:             target production used for enforced
+0026 %                                objetive limits (from minimum to maximum
+0027 %                                production)
+0028 %           - v_matrix:          fluxes for each target reaction predicted
+0029 %                                and each alpha.
+0030 %           - rxnTargets:        a list with all reactions with fluxes that
+0031 %                                change consistently as target production
+0032 %                                increases.
+0033 %                                Contains: ID, name, slope, gene rule, and
+0034 %                                equation
+0035 %           - transportTargets:  a list with all transport reactions with
+0036 %                                fluxes that change consistently as target
+0037 %                                production increases.
+0038 %                                Contains: ID, name, slope, gene rule, and
+0039 %                                equation
+0040 %           - geneTargets:       a list with all selected targets that
+0041 %                                increase production.
+0042 %                                Contains: gene, shortName, slope, action,
+0043 %                                and essentiality
 0044 %
 0045 % Usage:
-0046 %   FC = ecFSEOF(ecModel,targetRxn,csRxn,alphaLims,nSteps,filePath,filterG,modelAdapter)
+0046 %   fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,filterG,modelAdapter)
 0047 
-0048 if nargin < 8 || isempty(modelAdapter)
+0048 if nargin < 7 || isempty(modelAdapter)
 0049     modelAdapter = ModelAdapterManager.getDefault();
 0050     if isempty(modelAdapter)
 0051         error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
@@ -141,199 +141,270 @@ 

SOURCE CODE ^end 0054 params = modelAdapter.getParameters(); 0055 -0056 if nargin < 7 || isempty(filterG) -0057 filterG = true; +0056 if nargin < 5 || isempty(outputFile) +0057 outputFile = false; 0058 end 0059 0060 if nargin < 6 || isempty(filePath) 0061 filePath = fullfile(params.path,'output'); 0062 end 0063 -0064 if nargin < 5 || isempty(nSteps) +0064 if nargin < 4 || isempty(nSteps) 0065 nSteps = 16; 0066 end 0067 -0068 if nargin < 4 || isempty(alphaLims) -0069 alphaLims = [0.5 0.75]; -0070 end -0071 -0072 if numel(alphaLims) ~= 2 -0073 error('alphaLims parameter should be a vector of two values') -0074 end -0075 -0076 % Standardize grRules and rxnGeneMat in model -0077 [grRules,rxnGeneMat] = standardizeGrRules(ecModel,true); -0078 ecModel.grRules = grRules; -0079 ecModel.rxnGeneMat = rxnGeneMat; -0080 -0081 % Check carbon source uptake rate and LB -0082 ecModel = setParam(ecModel, 'obj', params.bioRxn, 1); -0083 sol = solveLP(ecModel, 1); -0084 csRxnIdx = strcmpi(ecModel.rxns,csRxn); +0068 % Get relevant rxn indexes +0069 prodTargetRxnIdx = getIndexes(model, prodTargetRxn,'rxns'); +0070 csRxnIdx = getIndexes(model, csRxn,'rxns'); +0071 bioRxnIdx = getIndexes(model, params.bioRxn,'rxns'); +0072 +0073 % Standardize grRules and rxnGeneMat in model +0074 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0075 model.grRules = grRules; +0076 model.rxnGeneMat = rxnGeneMat; +0077 +0078 % Check carbon source uptake rate and LB defined +0079 model = setParam(model, 'obj', params.bioRxn, 1); +0080 sol = solveLP(model); +0081 if model.lb(csRxnIdx) < sol.x(csRxnIdx) +0082 printOrange(['WARNING: Carbon source lower bound was set to ' num2str(model.lb(csRxnIdx)) ... +0083 ', but the uptake rate after model optimization is ' num2str(sol.x(csRxnIdx)) '.\n']) +0084 end 0085 -0086 if sol.x(csRxnIdx) < ecModel.lb(csRxnIdx) -0087 printOrange('WARNING: Carbon source LB and uptake rate are not equal.') -0088 end -0089 -0090 % run FSEOF analysis -0091 -0092 % Define alpha vector for suboptimal enforced objective values and -0093 % initialize fluxes and K_scores matrices -0094 alpha = alphaLims(1):((alphaLims(2)-alphaLims(1))/(nSteps-1)):alphaLims(2); -0095 v_matrix = zeros(length(ecModel.rxns),length(alpha)); -0096 k_matrix = zeros(length(ecModel.rxns),length(alpha)); -0097 enzRxns = contains(ecModel.rxns,'usage_prot_'); -0098 tolerance = 1e-4; -0099 -0100 % Simulate WT (100% growth) and forced (X% growth and the rest towards product): -0101 flux_WT = getFluxTarget(ecModel,params.bioRxn,targetRxn,1); +0086 % run FSEOF analysis +0087 +0088 % Find out the initial target production. +0089 iniTarget = sol.x(prodTargetRxnIdx); +0090 +0091 % Find out the maximum theoretical yield of target reaction. +0092 model = setParam(model,'obj',prodTargetRxn,1); +0093 sol = solveLP(model,1); +0094 % Set to 90%% based on https://doi.org/10.1128/AEM.00115-10 +0095 maxTarget = sol.x(prodTargetRxnIdx) * 0.9; +0096 +0097 % Define alpha vector for suboptimal enforced objective values between +0098 % minimal production and 90%% of the maximum theoretical yield, initialize +0099 % fluxes matriz +0100 alpha = iniTarget:((maxTarget-iniTarget)/(nSteps-1)):maxTarget; +0101 v_matrix = zeros(length(model.rxns),length(alpha)); 0102 -0103 % Set values which are under solver tolerance -0104 flux_WT(flux_WT < 0 & ~enzRxns) = 0; -0105 flux_WT(flux_WT > 0 & enzRxns) = 0; -0106 -0107 progressbar('Flux Scanning with Enforced Objective Function') -0108 for i = 1:length(alpha) -0109 flux_MAX = getFluxTarget(ecModel,params.bioRxn,targetRxn,alpha(i)); -0110 % Set values which are under solver tolerance -0111 flux_MAX(flux_MAX < 0 & ~enzRxns) = 0; -0112 flux_MAX(flux_MAX > 0 & enzRxns) = 0; -0113 v_matrix(:,i) = flux_MAX; -0114 k_matrix(:,i) = flux_MAX./flux_WT; -0115 progressbar(i/length(alpha)) -0116 end -0117 progressbar(1) % Make sure it closes -0118 -0119 % take out rxns with no grRule and standard gene -0120 withGR = ~cellfun(@isempty,ecModel.grRules); -0121 stdPos = contains(ecModel.grRules,'standard'); -0122 withGR(stdPos) = 0; -0123 rxns = ecModel.rxns(withGR); -0124 v_matrix = v_matrix(withGR,:); -0125 k_matrix = k_matrix(withGR,:); -0126 rxnGeneM = ecModel.rxnGeneMat(withGR,:); -0127 -0128 % Filter out rxns that are always zero -> k=0/0=NaN: -0129 non_nan = sum(~isnan(k_matrix),2) > 0; -0130 rxns = rxns(non_nan,:); -0131 v_matrix = v_matrix(non_nan,:); -0132 k_matrix = k_matrix(non_nan,:); -0133 rxnGeneM = rxnGeneM(non_nan,:); +0103 % Enforce objective flux iteratively +0104 progressbar('Flux Scanning with Enforced Objective Function') +0105 for i = 1:nSteps +0106 % Enforce the objetive flux of product formation +0107 model = setParam(model,'eq',prodTargetRxnIdx,alpha(i)); +0108 +0109 % Restore minimum biomass lb to zero and set it as objective +0110 model.lb(bioRxnIdx) = 0; +0111 model = setParam(model,'obj',params.bioRxn,1); +0112 sol = solveLP(model,1); +0113 +0114 % Store flux distribution +0115 v_matrix(:,i) = sol.x; +0116 +0117 progressbar(i/nSteps) +0118 end +0119 progressbar(1) % Make sure it closes +0120 +0121 % Take out rxns with no grRule and standard gene +0122 withGR = ~cellfun(@isempty,model.grRules); +0123 stdIdx = contains(model.grRules,'standard'); +0124 withGR(stdIdx) = 0; +0125 target_rxns = model.rxns(withGR); +0126 v_matrix = v_matrix(withGR,:); +0127 rxnGeneM = model.rxnGeneMat(withGR,:); +0128 +0129 % Filter out rxns that are always zero +0130 zero_flux = ~all(abs(v_matrix(:,1:nSteps))==0,2); +0131 target_rxns = target_rxns(zero_flux,:); +0132 v_matrix = v_matrix(zero_flux,:); +0133 rxnGeneM = rxnGeneM(zero_flux,:); 0134 -0135 % Replace remaining NaNs with 1s: -0136 k_matrix(isnan(k_matrix)) = 1; -0137 -0138 % Replace any Inf value with 1000 (maximum value is ~700): -0139 k_matrix(k_matrix > 1000) = 1000; -0140 % k_matrix(k_matrix < -1000) = 2; -0141 -0142 % Filter out values that are inconsistent at different alphas: -0143 always_down = sum(k_matrix <= 1,2) == length(alpha); -0144 always_up = sum(k_matrix >= 1,2) == length(alpha); -0145 -0146 % Identify those reactions with mixed patterns -0147 incons_rxns = always_down + always_up == 0; -0148 -0149 % Identify genes that are linked to "inconsistent rxns" -0150 incons_genes = sum(rxnGeneM(incons_rxns,:),1) > 0; -0151 -0152 % Finally, inconsistent reactions are those that are not conected -0153 % to "inconsistent genes" from the original "inconsistent rxns" set -0154 incons_rxns = sum(rxnGeneM(:,incons_genes),2) > 0; -0155 -0156 % Keep results for the consistent rxns exclusively -0157 rxns = rxns(~incons_rxns,:); -0158 v_matrix = v_matrix(~incons_rxns,:); -0159 k_matrix = k_matrix(~incons_rxns,:); -0160 rxnGeneM = rxnGeneM(~incons_rxns,:); -0161 -0162 % Get median k-score across steps and order from highest to lowest -0163 k_rxns = mean(k_matrix,2); -0164 [~,order] = sort(k_rxns,'descend'); -0165 rxns = rxns(order,:); -0166 v_matrix = v_matrix(order,:); -0167 k_matrix = k_matrix(order,:); -0168 rxnGeneM = rxnGeneM(order,:); -0169 k_rxns = k_rxns(order,:); -0170 -0171 % Create list of remaining genes and filter out any inconsistent score: -0172 % Just those genes that are connected to the remaining rxns are -0173 genes = ecModel.genes(sum(rxnGeneM,1) > 0); -0174 k_genes = zeros(size(genes)); -0175 cons_genes = false(size(genes)); -0176 rxnGeneM = rxnGeneM(:,sum(rxnGeneM,1) > 0); -0177 -0178 for i = 1:length(genes) -0179 % Extract all the K_scores (from rxns across alphas) conected to -0180 % each remaining gene -0181 k_set = k_rxns(rxnGeneM(:,i) > 0); -0182 % Check the kind of control that gene i-th exerts over its reactions -0183 always_down = sum(k_set <= (1-tolerance)) == length(k_set); -0184 always_up = sum(k_set >= (1+tolerance)) == length(k_set); -0185 % Evaluate if gene is always exerting either a positive or negative -0186 % control -0187 cons_genes(i) = always_down + always_up == 1; -0188 k_genes(i) = mean(k_set); -0189 end -0190 -0191 % Keep "consistent genes" -0192 genes = genes(cons_genes); -0193 k_genes = k_genes(cons_genes); -0194 rxnGeneM = rxnGeneM(:,cons_genes); +0135 % Identify those rxns that always increase or decrease, and calculate the +0136 % slope as the difference in the flux when enforce objetive target +0137 % production is set to 90%% of the maximum teorethical yield +0138 % << v_matrix(i,nSteps-1) >> and the flux when the enforce objetive target +0139 % production is set to the minimum << v_matrix(i,1) >> for the reaction i, +0140 % divided by maxTarget-maxTarget/nSteps-1. +0141 slope_rxns = zeros(size(target_rxns)); +0142 targets = logical(size(target_rxns)); +0143 target_type = cell(size(target_rxns)); +0144 for i = 1:length(target_rxns) +0145 if issorted(abs(v_matrix(i,1:nSteps)),'strictascend') +0146 % Those reactions that always increase while enforcing target +0147 % production are suggested for Over Expression +0148 targets(i) = true; +0149 slope_rxns(i) = abs(v_matrix(i,nSteps-1)-v_matrix(i,1))/abs(maxTarget-maxTarget/nSteps-1); +0150 target_type(i) = {'OE'}; +0151 elseif issorted(abs(v_matrix(i,1:nSteps)),'strictdescend') +0152 % Those reactions that always decrease while enforcing target +0153 % production are suggested for KnockDown or KnockOut. KO are those +0154 % reactions which have zero flux when enforcing target production +0155 % to 90%% of the maximum theoretical yield. +0156 targets(i) = true; +0157 slope_rxns(i) = abs(v_matrix(i,nSteps-1)-v_matrix(i,1))/abs(maxTarget-maxTarget/nSteps-1); +0158 if v_matrix(i,nSteps) == 0 +0159 target_type(i) = {'KO'}; +0160 else +0161 target_type(i) = {'KD'}; +0162 end +0163 end +0164 end +0165 +0166 % Only keep those reaction that shows an increase or decrease pattern. +0167 target_rxns = target_rxns(targets); +0168 v_matrix = v_matrix(targets,:); +0169 rxnGeneM = rxnGeneM(targets,:); +0170 slope_rxns = slope_rxns(targets); +0171 target_type = target_type(targets); +0172 +0173 % Order from highest to lowest slope +0174 [~,order] = sort(slope_rxns,'descend'); +0175 target_rxns = target_rxns(order); +0176 v_matrix = v_matrix(order,:); +0177 rxnGeneM = rxnGeneM(order,:); +0178 slope_rxns = slope_rxns(order); +0179 target_type = target_type(order); +0180 +0181 % Filter out reactions with slope < quantile +0182 non_zero_slope = slope_rxns > quantile(slope_rxns,0.75); +0183 target_rxns = target_rxns(non_zero_slope); +0184 v_matrix = v_matrix(non_zero_slope,:); +0185 rxnGeneM = rxnGeneM(non_zero_slope,:); +0186 slope_rxns = slope_rxns(non_zero_slope); +0187 target_type = target_type(non_zero_slope); +0188 +0189 % Create gene list of those connected to the remaining rxns +0190 genes = model.genes(sum(rxnGeneM,1) > 0); +0191 slope_genes = zeros(size(genes)); +0192 rxnGeneM = rxnGeneM(:,sum(rxnGeneM,1) > 0); +0193 target_type_genes = cell(size(genes)); +0194 essentiality = cell(size(genes)); 0195 -0196 if filterG -0197 % Filter any value between mean(alpha) and 1: -0198 unchanged = (k_genes >= mean(alpha) - tolerance) + (k_genes <= 1 + tolerance) == 2; -0199 genes = genes(~unchanged); -0200 k_genes = k_genes(~unchanged); -0201 rxnGeneM = rxnGeneM(:,~unchanged); -0202 % Update results for rxns-related fields (remove remaining reactions -0203 % without any associated gene in rxnGeneM) -0204 toKeep = (sum(rxnGeneM,2) > 0); -0205 rxns = rxns(toKeep,:); -0206 v_matrix = v_matrix(toKeep,:); -0207 k_matrix = k_matrix(toKeep,:); -0208 k_rxns = k_rxns(toKeep,:); -0209 end -0210 -0211 % Order genes from highest to lowest k: -0212 [~,order] = sort(k_genes,'descend'); -0213 genes = genes(order,:); -0214 k_genes = k_genes(order,:); -0215 -0216 % Create output (exclude enzyme usage reactions): -0217 toKeep = ~startsWith(rxns,'usage_prot_'); -0218 rxnIdx = getIndexes(ecModel,rxns(toKeep),'rxns'); -0219 geneIdx = cellfun(@(x) find(strcmpi(ecModel.genes,x)),genes); -0220 FC.flux_WT = flux_WT; -0221 FC.alpha = alpha; -0222 FC.v_matrix = v_matrix(toKeep,:); -0223 FC.k_matrix = k_matrix(toKeep,:); -0224 FC.rxnsTable(:,1) = ecModel.rxns(rxnIdx); -0225 FC.rxnsTable(:,2) = ecModel.rxnNames(rxnIdx); -0226 FC.rxnsTable(:,3) = num2cell(k_rxns(toKeep)); -0227 FC.rxnsTable(:,4) = ecModel.grRules(rxnIdx); -0228 FC.rxnsTable(:,5) = constructEquations(ecModel,rxnIdx); -0229 FC.geneTable(:,1) = genes; -0230 FC.geneTable(:,2) = ecModel.geneShortNames(geneIdx); -0231 FC.geneTable(:,3) = num2cell(k_genes); -0232 -0233 writetable(cell2table(FC.geneTable, ... -0234 'VariableNames', {'gene_IDs' 'gene_names' 'K_score'}), ... -0235 fullfile(filePath, 'ecFSEOF_genes.tsv'), ... -0236 'FileType', 'text', ... -0237 'Delimiter', '\t', ... -0238 'QuoteStrings', false); -0239 -0240 writetable(cell2table(FC.rxnsTable, ... -0241 'VariableNames', {'rxn_IDs' 'rxnNames' 'K_score' 'grRules' 'rxn_formula'}), ... -0242 fullfile(filePath, 'ecFSEOF_rxns.tsv'), ... -0243 'FileType', 'text', ... -0244 'Delimiter', '\t', ... -0245 'QuoteStrings', false); -0246 -0247 disp(['ecFSEOF results stored at: ' newline filePath]); -0248 end

+0196 % Validate for gene essentiality +0197 progressbar('Checking for gene essentiality') +0198 for i = 1:numel(genes) +0199 +0200 % Block protein usage to KO al the reactions associated to it. +0201 usage_rxn_idx = strcmpi(model.ec.genes, genes{i}); +0202 usage_rxn = strcat('usage_prot_', model.ec.enzymes(usage_rxn_idx)); +0203 tempModel = setParam(model, 'eq', usage_rxn, 0); +0204 solKO = solveLP(tempModel); +0205 % Check if no feasible solution was found +0206 if solKO.stat == -1 || solKO.x(bioRxnIdx) < 1e-8 +0207 essentiality(i) = {'essential'}; +0208 end +0209 +0210 % Since a gene can be involved in multiple reactions, multiple +0211 % engineering manipulations can be suggested for the same gene +0212 % e.g. (OE and KD). So, reconcilie them, and report only one. +0213 +0214 % Get reaction index, in targets set, for each gene +0215 rxns_for_gene = find(rxnGeneM(:,i) > 0); +0216 actions = unique(target_type(rxns_for_gene)); +0217 if numel(actions) > 1 +0218 % If any OE is suggested give the highest priority over KD or KO +0219 if any(startsWith(actions, 'OE')) +0220 target_type_genes(i) = {'OE'}; +0221 % Otherwise change to KO if not essential +0222 elseif ~any(startsWith(actions, 'OE')) && ~isequal(essentiality(i),{'essential'}) +0223 target_type_genes(i) = {'KO'}; +0224 else +0225 target_type_genes(i) = {'KD'}; +0226 end +0227 else +0228 target_type_genes(i) = actions; +0229 end +0230 % Extract all the slope (from rxns across alphas) conected to +0231 % each remaining gene +0232 slope_set = slope_rxns(rxns_for_gene); +0233 slope_genes(i) = mean(slope_set); +0234 +0235 progressbar(i/numel(genes)) +0236 end +0237 progressbar(1) % Make sure it closes +0238 +0239 % Order genes from highest to lowest slope: +0240 [~,order] = sort(slope_genes,'descend'); +0241 genes = genes(order,:); +0242 slope_genes = slope_genes(order,:); +0243 target_type_genes = target_type_genes(order,:); +0244 essentiality = essentiality(order,:); +0245 +0246 % Create output: +0247 +0248 % Report values used to enforce objective target production +0249 fseof.alpha = alpha; +0250 +0251 % Exclude enzyme usage reactions +0252 toKeep = ~startsWith(target_rxns,'usage_prot_'); +0253 rxnIdx = getIndexes(model,target_rxns(toKeep),'rxns'); +0254 +0255 % Report the fluxes at each alpha +0256 fseof.v_matrix = array2table(v_matrix(toKeep,:), ... +0257 'VariableNames', string(alpha), 'RowNames', model.rxns(rxnIdx)); +0258 +0259 % Report reaction targets +0260 fseof.rxnTargets(:,1) = model.rxns(rxnIdx); +0261 fseof.rxnTargets(:,2) = model.rxnNames(rxnIdx); +0262 fseof.rxnTargets(:,3) = num2cell(slope_rxns(toKeep)); +0263 fseof.rxnTargets(:,4) = model.grRules(rxnIdx); +0264 fseof.rxnTargets(:,5) = constructEquations(model,rxnIdx); +0265 +0266 % Identify transport reactions; defined as reactions involving (at least) +0267 % one metabolite name in more than one compartment. +0268 transportRxns = false(numel(rxnIdx),1); +0269 for i = 1:numel(rxnIdx) +0270 % Get the involved metabolites in each reaction +0271 mets = model.metNames(model.S(:,rxnIdx(i))~=0); +0272 % Remove enzyme metabolite +0273 mets(startsWith(mets,'prot_')) = []; +0274 % Validate if it is a transport reaciton +0275 transportRxns(i) = numel(mets) ~= numel(unique(mets)); +0276 end +0277 +0278 % Create separate table for transport reactions +0279 fseof.transportTargets = fseof.rxnTargets(transportRxns,:); +0280 fseof.rxnTargets(transportRxns,:) = []; +0281 +0282 % Report gene targets +0283 geneIdx = cellfun(@(x) find(strcmpi(model.genes,x)),genes); +0284 fseof.geneTargets(:,1) = genes; +0285 fseof.geneTargets(:,2) = model.geneShortNames(geneIdx); +0286 fseof.geneTargets(:,3) = num2cell(slope_genes); +0287 fseof.geneTargets(:,4) = target_type_genes; +0288 fseof.geneTargets(:,5) = essentiality; +0289 +0290 % Save results in a file if defined +0291 if outputFile +0292 % Write file with gene targets +0293 writetable(cell2table(fseof.geneTargets, ... +0294 'VariableNames', {'gene_IDs' 'gene_names' 'slope' 'action' 'essentiality'}), ... +0295 fullfile(filePath, 'ecFSEOF_genes.tsv'), ... +0296 'FileType', 'text', ... +0297 'Delimiter', '\t', ... +0298 'QuoteStrings', false); +0299 +0300 % Write file with rxn targets +0301 writetable(cell2table(fseof.rxnTargets, ... +0302 'VariableNames', {'rxn_IDs' 'rxnNames' 'slope' 'grRules' 'rxn_formula'}), ... +0303 fullfile(filePath, 'ecFSEOF_rxns.tsv'), ... +0304 'FileType', 'text', ... +0305 'Delimiter', '\t', ... +0306 'QuoteStrings', false); +0307 +0308 % Write file with transport targets +0309 writetable(cell2table(fseof.transportTargets, ... +0310 'VariableNames', {'rxn_IDs' 'rxnNames' 'slope' 'grRules' 'rxn_formula'}), ... +0311 fullfile(filePath, 'ecFSEOF_transporter.tsv'), ... +0312 'FileType', 'text', ... +0313 'Delimiter', '\t', ... +0314 'QuoteStrings', false); +0315 +0316 disp(['ecFSEOF results stored at: ' newline filePath]); +0317 end +0318 +0319 end

Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/src/geckomat/utilities/getFluxTarget.html b/doc/src/geckomat/utilities/getFluxTarget.html deleted file mode 100644 index f790eec11..000000000 --- a/doc/src/geckomat/utilities/getFluxTarget.html +++ /dev/null @@ -1,133 +0,0 @@ - - - - Description of getFluxTarget - - - - - - - - - -
Home > src > geckomat > utilities > getFluxTarget.m
- - - -

getFluxTarget -

- -

PURPOSE ^

-
getFluxTarget
- -

SYNOPSIS ^

-
function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)
- -

DESCRIPTION ^

-
 getFluxTarget
-
-   Function that performs a series of LP optimizations on an ecModel,
-   by first maximizing biomass, then fixing a suboptimal value and 
-   proceeding to protein pool minimization, last a minimization and
-   maximization of a given production target reaction is performed
-   subject to the optimal production level.
-
- Input:
-   model           an ecModel in GECKO 3 format (with ecModel.ec structure).
-   bioRxn          rxn ID for the biomass reaction.
-   targetRxn       rxn ID for the production target reaction, a exchange
-                   reaction is recommended.
-   alpha           scalling factor for desired suboptimal growth.
-                   (Optional, default 0.95)
-
- Output:
-   minFlux         vector of minimum flux rates at minimum target production,
-                   corresponding to model.rxns
-   maxFlux         vector of maximum flux rates at maximum target production,
-                   corresponding to model.rxns
-
- Usage:
-   [maxFlux, minFlux] = getFluxTarget(ecModel,bioRxn,targetRxn,alpha)
- - -

CROSS-REFERENCE INFORMATION ^

-This function calls: -
    -
-This function is called by: - - - - - -

SOURCE CODE ^

-
0001 function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)
-0002 % getFluxTarget
-0003 %
-0004 %   Function that performs a series of LP optimizations on an ecModel,
-0005 %   by first maximizing biomass, then fixing a suboptimal value and
-0006 %   proceeding to protein pool minimization, last a minimization and
-0007 %   maximization of a given production target reaction is performed
-0008 %   subject to the optimal production level.
-0009 %
-0010 % Input:
-0011 %   model           an ecModel in GECKO 3 format (with ecModel.ec structure).
-0012 %   bioRxn          rxn ID for the biomass reaction.
-0013 %   targetRxn       rxn ID for the production target reaction, a exchange
-0014 %                   reaction is recommended.
-0015 %   alpha           scalling factor for desired suboptimal growth.
-0016 %                   (Optional, default 0.95)
-0017 %
-0018 % Output:
-0019 %   minFlux         vector of minimum flux rates at minimum target production,
-0020 %                   corresponding to model.rxns
-0021 %   maxFlux         vector of maximum flux rates at maximum target production,
-0022 %                   corresponding to model.rxns
-0023 %
-0024 % Usage:
-0025 %   [maxFlux, minFlux] = getFluxTarget(ecModel,bioRxn,targetRxn,alpha)
-0026 
-0027 if nargin < 4 || isempty(alpha)
-0028     alpha = 0.95;
-0029 end
-0030 
-0031 % Fix a suboptimal alpha if equal to 1
-0032 if alpha == 1
-0033     alpha = 0.99;
-0034 end
-0035 
-0036 % Get relevant rxn indexes
-0037 bioRxnIdx    = getIndexes(model, bioRxn,'rxns');
-0038 poolIdx      = strcmpi(model.rxns, 'prot_pool_exchange');
-0039 
-0040 % Maximize growth and fix carbon source and suboptimal growth
-0041 model = setParam(model, 'obj', bioRxn, 1);
-0042 sol   = solveLP(model);
-0043 model = setParam(model, 'lb', bioRxn, sol.x(bioRxnIdx) * alpha);
-0044 
-0045 % Minimize prot_pool_exchange and fix
-0046 model = setParam(model, 'obj', 'prot_pool_exchange', 1);
-0047 sol   = solveLP(model);
-0048 model = setParam(model, 'lb', 'prot_pool_exchange', sol.x(poolIdx) * 1.01);
-0049 
-0050 % If minimum fluxes are required get them
-0051 minFlux = [];
-0052 if nargout == 2
-0053     % Minimize target
-0054     model = setParam(model, 'obj', targetRxn, -1);
-0055     minSol = solveLP(model);
-0056     minFlux = minSol.x;
-0057 end
-0058 
-0059 % Maximize target
-0060 model = setParam(model, 'obj', targetRxn, 1);
-0061 maxSol = solveLP(model);
-0062 maxFlux = maxSol.x;
-0063 end
-
Generated by m2html © 2005
- - \ No newline at end of file diff --git a/doc/src/geckomat/utilities/index.html b/doc/src/geckomat/utilities/index.html index 49d98ee11..fa531e367 100644 --- a/doc/src/geckomat/utilities/index.html +++ b/doc/src/geckomat/utilities/index.html @@ -19,7 +19,7 @@

Index for src\geckomat\utilities

Matlab files in this directory:

-
 addNewRxnsToECaddNewRxnsToEC
 ecFSEOFecFSEOF
 ecFVAecFVA
 enzymeUsageenzymeUsage
 findGECKOrootfindGECKOroot
 getFluxTargetgetFluxTarget
 getSubsetEcModelgetSubset_ecModel
 loadConventionalGEMloadConventionalGEM
 loadEcModelloadEcModel
 loadFluxDataloadFluxData
 loadProtDataloadProtData
 mapRxnsToConvmapRxnsToConv
 plotEcFVAplotEcFVA
 reportEnzymeUsagereportEnzymeUsage
 saveEcModelsaveECmodel
 startGECKOprojectstartGECKOproject
 updateGECKOdocupdateGeckoDoc
+ addNewRxnsToECaddNewRxnsToEC  ecFSEOFecFSEOF  ecFVAecFVA  enzymeUsageenzymeUsage  findGECKOrootfindGECKOroot  getSubsetEcModelgetSubset_ecModel  loadConventionalGEMloadConventionalGEM  loadEcModelloadEcModel  loadFluxDataloadFluxData  loadProtDataloadProtData  mapRxnsToConvmapRxnsToConv  plotEcFVAplotEcFVA  reportEnzymeUsagereportEnzymeUsage  saveEcModelsaveECmodel  startGECKOprojectstartGECKOproject  updateGECKOdocupdateGeckoDoc diff --git a/doc/src/geckomat/utilities/startGECKOproject.html b/doc/src/geckomat/utilities/startGECKOproject.html index 2926ba10b..f51e357ec 100644 --- a/doc/src/geckomat/utilities/startGECKOproject.html +++ b/doc/src/geckomat/utilities/startGECKOproject.html @@ -113,7 +113,8 @@

SOURCE CODE ^else 0058 printOrange('WARNING: A project with the same name exits at the same location. The project was not created.\n') 0059 end -0060 end

+0060 cd(fullPath) +0061 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/src/geckomat/change_model/makeEcModel.m b/src/geckomat/change_model/makeEcModel.m index 1bd53fdab..5fabb4bc1 100644 --- a/src/geckomat/change_model/makeEcModel.m +++ b/src/geckomat/change_model/makeEcModel.m @@ -326,7 +326,7 @@ %9: Add proteins as pseudometabolites if ~geckoLight - [proteinMets.mets, uniprotSortId] = sort(ec.enzymes); + [proteinMets.mets, uniprotSortId] = unique(ec.enzymes); proteinMets.mets = strcat('prot_',proteinMets.mets); proteinMets.metNames = proteinMets.mets; proteinMets.compartments = compartmentID; diff --git a/src/geckomat/utilities/startGECKOproject.m b/src/geckomat/utilities/startGECKOproject.m index 4b85b5fc0..f6f2c2c88 100644 --- a/src/geckomat/utilities/startGECKOproject.m +++ b/src/geckomat/utilities/startGECKOproject.m @@ -57,4 +57,5 @@ function startGECKOproject(name, path) else printOrange('WARNING: A project with the same name exits at the same location. The project was not created.\n') end +cd(fullPath) end diff --git a/tutorials/light_ecModel/HumanGEMAdapter.m b/tutorials/light_ecModel/HumanGEMAdapter.m index a7e393237..4a9ae97ba 100644 --- a/tutorials/light_ecModel/HumanGEMAdapter.m +++ b/tutorials/light_ecModel/HumanGEMAdapter.m @@ -6,7 +6,7 @@ % The model distributed with the light_ecModel tutorial is Human-GEM % is version 1.15.0, available from - % https://github.com/SysBioChalmers/Human-GEM/releases/tag/v1.3.0 + % https://github.com/SysBioChalmers/Human-GEM/releases/tag/v1.15.0 % In addition, the following lines were run to reduce its size % before storing it in this GECKO tutorial: % ihuman = simplifyModel(ihuman,false,false,true,true);