ecFSEOF

 ecFSEOF
   Function that runs Flux-Scanning with Enforced Objective Function (FSEOF)
   for a specified production target.

 Input:
   model          an ecModel in GECKO 3 format (with ecModel.ec structure).
   prodTargetRxn  rxn ID for the production target reaction, a exchange
                  reaction is recommended.
   csRxn          rxn ID for the main carbon source uptake reaction.
   nSteps         number of steps for suboptimal objective in FSEOF.
                  (Optional, default 16)
   outputFile     bolean option to save results in a file. (Optional,
                  default false)
   filePath       file path for results output. It will store two files:
                  - at the genes level, ecFSEOF_genes.tsv
                  - at the reactions level, ecFSEOF_rxns.tsv
                  (Optional, default in the 'output' sub-folder taken from
                  modelAdapter, e.g. output/ecFSEOF_rxns.tsv)
   modelAdapter   a loaded model adapter. (Optional, will otherwise use
                  the default model adapter)

 Output:
   fseof   an structure with all results. Contains the following fields:
           - alpha:             target production used for enforced 
                                objetive limits (from minimum to maximum 
                                production)
           - v_matrix:          fluxes for each target reaction predicted
                                and each alpha.
           - rxnTargets:        a list with all reactions with fluxes that
                                change consistently as target production
                                increases.
                                Contains: ID, name, slope, gene rule, and
                                equation
           - transportTargets:  a list with all transport reactions with
                                fluxes that change consistently as target
                                production increases.
                                Contains: ID, name, slope, gene rule, and
                                equation
           - geneTargets:       a list with all selected targets that
                                increase production.
                                Contains: gene, shortName, slope, action,
                                and essentiality

 Usage:
   fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,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 %   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 %   fseof = ecFSEOF(model,prodTargetRxn,csRxn,nSteps,outputFile,filePath,filterG,modelAdapter)
0047 
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.')
0052     end
0053 end
0054 params = modelAdapter.getParameters();
0055 
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 < 4 || isempty(nSteps)
0065     nSteps = 16;
0066 end
0067 
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 % 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 % 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 % 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 % 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