addNewRxnsToEC

 addNewRxnsToEC
   Add new reaction to an enzyme-constrained model. This function is 
   useful to simulate metabolic manipulations done to an organism such as
   integration of new genes that add a new reaction/pathway.

 Input:
   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
   newRxns         structure with the new reaction information as follow:
                   rxns        cell array with unique strings that
                               identifies each reaction
                   rxnNames    cell array with the names of each reaction
                   equations   cell array with equation strings. Decimal
                               coefficients are expressed as "1.2".
                               Reversibility is indicated by "<=>" or
                               "=>".
                   grRules     cell array with the gene-reaction
                               relationship for each reaction. For example
                               "(A and B) or (C)" means that the reaction
                               could be catalyzed by a complex between
                               A & B or by C on its own. All the genes
                               have to be present in model.genes. Add
                               genes with addGenesRaven before calling
                               this function if needed (opt, default '')
   newEnzymes      structure with the new enzymes information as follow:
                   enzymes     cell array with uniprot id
                   genes       cell array with the respective gene
                   mw          cell array with the MW
   modelAdapter    a loaded model adapter (Optional, will otherwise use the
                   default model adapter).

 Output:
   model           ecModel whit new reactions
   rxnsAdded       cell array with the reactions added

 Notes:
     (i) Write equations as follows:
               'xylitol[c] + NAD[c] => D-xylulose[c] + NADH[c] + H+[c]'
         Note that the reversibility defined in the equation will be used
         to split to construct irreversible reactions (add _REV) and the
         rules defined will be used to expand the model (add _EXP_n).
        
     (2) After add the new reactions, setKcatForReactions or
         applyCustomKcats should be run to add the kcat to the reactions,
         and subsequently applyKcatConstraints.

 Usage:
   [model, rxnsAdded, enzAdded] = addNewRxnsToEC(model, newRxns, newEnzymes, modelAdapter);

0001 function [model, rxnsAdded, enzAdded] = addNewRxnsToEC(model, newRxns, newEnzymes, modelAdapter)
0002 % addNewRxnsToEC
0003 %   Add new reaction to an enzyme-constrained model. This function is
0004 %   useful to simulate metabolic manipulations done to an organism such as
0005 %   integration of new genes that add a new reaction/pathway.
0006 %
0007 % Input:
0008 %   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
0009 %   newRxns         structure with the new reaction information as follow:
0010 %                   rxns        cell array with unique strings that
0011 %                               identifies each reaction
0012 %                   rxnNames    cell array with the names of each reaction
0013 %                   equations   cell array with equation strings. Decimal
0014 %                               coefficients are expressed as "1.2".
0015 %                               Reversibility is indicated by "<=>" or
0016 %                               "=>".
0017 %                   grRules     cell array with the gene-reaction
0018 %                               relationship for each reaction. For example
0019 %                               "(A and B) or (C)" means that the reaction
0020 %                               could be catalyzed by a complex between
0021 %                               A & B or by C on its own. All the genes
0022 %                               have to be present in model.genes. Add
0023 %                               genes with addGenesRaven before calling
0024 %                               this function if needed (opt, default '')
0025 %   newEnzymes      structure with the new enzymes information as follow:
0026 %                   enzymes     cell array with uniprot id
0027 %                   genes       cell array with the respective gene
0028 %                   mw          cell array with the MW
0029 %   modelAdapter    a loaded model adapter (Optional, will otherwise use the
0030 %                   default model adapter).
0031 %
0032 % Output:
0033 %   model           ecModel whit new reactions
0034 %   rxnsAdded       cell array with the reactions added
0035 %
0036 % Notes:
0037 %     (i) Write equations as follows:
0038 %               'xylitol[c] + NAD[c] => D-xylulose[c] + NADH[c] + H+[c]'
0039 %         Note that the reversibility defined in the equation will be used
0040 %         to split to construct irreversible reactions (add _REV) and the
0041 %         rules defined will be used to expand the model (add _EXP_n).
0042 %
0043 %     (2) After add the new reactions, setKcatForReactions or
0044 %         applyCustomKcats should be run to add the kcat to the reactions,
0045 %         and subsequently applyKcatConstraints.
0046 %
0047 % Usage:
0048 %   [model, rxnsAdded, enzAdded] = addNewRxnsToEC(model, newRxns, newEnzymes, modelAdapter);
0049 
0050 if nargin < 4 || isempty(modelAdapter)
0051     modelAdapter = ModelAdapterManager.getDefault();
0052     if isempty(modelAdapter)
0053         error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
0054     end
0055 end
0056 params = modelAdapter.getParameters();
0057 
0058 if nargin < 3 || isempty(newEnzymes)
0059     newEnzymes = [];
0060 end
0061 
0062 % Validate if GECKO version 3
0063 if ~isfield(model, 'ec')
0064     error('Input model is not an ecModel in GECKO 3 format (with ecModel.ec structure)')
0065 end
0066 
0067 if model.ec.geckoLight
0068     error('addNewRxnsToEC only works on non-light ecModels')
0069 else
0070     % Check if the new proteins are not already present in the model
0071     if ~isempty(newEnzymes)
0072         [toRemove,~] = ismember(newEnzymes.enzymes, model.ec.enzymes);
0073 
0074         % Remove from the list those enzymes that already present
0075         if any(toRemove)
0076             printOrange(['WARNING: Enzymes ' strjoin(newEnzymes.enzymes(toRemove),', ') ' are already present in the model and will not be added.\n']);
0077             newEnzymes.enzymes(toRemove) = [];
0078             newEnzymes.genes(toRemove) = [];
0079             newEnzymes.mw(toRemove) = [];
0080         end
0081 
0082         if ~isempty(newEnzymes.enzymes)
0083             
0084             % 1. Add new genes
0085             genesToAdd.genes = newEnzymes.genes;
0086             genesToAdd.geneShortNames = genesToAdd.genes;
0087             model = addGenesRaven(model, genesToAdd);
0088 
0089             % 2. Add new Enzyme-usage pseudometabolite
0090 
0091             % Validate compartment to add protein pseudoreactions
0092             compartmentID = strcmp(model.compNames,params.enzyme_comp);
0093             if ~any(compartmentID)
0094                 error(['Compartment ' params.enzyme_comp ' (specified in params.enzyme_comp) '...
0095                     'cannot be found in model.compNames'])
0096             end
0097             compartmentID = model.comps(compartmentID);
0098 
0099             proteinMets.mets         = strcat('prot_', newEnzymes.enzymes);
0100             proteinMets.metNames     = proteinMets.mets;
0101             proteinMets.compartments = compartmentID;
0102             if isfield(model,'metMiriams')
0103                 proteinMets.metMiriams = repmat({struct('name',{{'sbo'}},'value',{{'SBO:0000252'}})},numel(proteinMets.mets),1);
0104             end
0105             if isfield(model,'metCharges')
0106                 proteinMets.metCharges = zeros(numel(proteinMets.mets),1);
0107             end
0108             proteinMets.metNotes = repmat({'Enzyme-usage pseudometabolite'},numel(proteinMets.mets),1);
0109             model = addMets(model,proteinMets);
0110         end
0111         enzAdded = newEnzymes.enzymes(:);  
0112     end
0113 
0114     % 3. Add new reactions
0115 
0116     % Before add the reactions validate that all the genes in grRules are
0117     % present in the model.
0118     genesInRules = strjoin(convertCharArray(newRxns.grRules));
0119     genesInRules = regexp(genesInRules,' |)|(|and|or','split'); % Remove all grRule punctuation
0120     genesInRules = genesInRules(~cellfun(@isempty,genesInRules));  % Remove spaces and empty genes
0121     genesInRules = setdiff(unique(genesInRules),model.genes); % Only keep new genes
0122 
0123     % If there is any gene which is not present, as to be pass as
0124     % newEnzymes input
0125     if ~isempty(genesInRules)
0126         errorText = 'The following genes from grRules are not present in the model, add them as newEnzymes input:';
0127         dispEM(errorText,true,genesInRules,false)
0128     end
0129 
0130     % Get reversible reactions to split
0131     toIrrev = find(cellfun(@(x) contains(x,'<=>'),newRxns.equations));
0132     if ~isempty(toIrrev)
0133         for i=1:numel(toIrrev)
0134             idx = toIrrev(i);
0135             newEquation = split(newRxns.equations(idx),'<=>');
0136             % Remove directionality
0137             newRxns.equations(idx) = erase(newRxns.equations(idx), '<');
0138 
0139             % Add a new reaction to the list to be added
0140             newRxns.rxns{end+1} = [newRxns.rxns{idx} '_REV'];
0141             newRxns.rxnNames{end+1} = [newRxns.rxnNames{idx} ' (reversible)'];
0142             newRxns.equations{end+1} = [strtrim(newEquation{2, 1}) ' => ' strtrim(newEquation{1, 1})];
0143             newRxns.grRules{end+1} = newRxns.grRules{idx};
0144         end
0145     end
0146 
0147     % Get reactions to separate isozymes
0148     toSplit = find(cellfun(@(x) contains(x,'or'),newRxns.grRules));
0149     if ~isempty(toSplit)
0150         for i=1:numel(toSplit)
0151             idx = toSplit(i);
0152             newRules = split(newRxns.grRules(idx),'or');
0153             newRules = strtrim(erase(newRules, {'(' ')'}));
0154 
0155             for j=1:numel(newRules)
0156                 newRxns.rxns{end+1} = [newRxns.rxns{idx} '_EXP_' num2str(j)];
0157                 newRxns.rxnNames(end+1) = newRxns.rxnNames(idx);
0158                 newRxns.equations(end+1) = newRxns.equations(idx);
0159                 newRxns.grRules(end+1) = newRules(j);
0160             end
0161 
0162             % Remove the original one
0163             newRxns.rxns(idx) = [];
0164             newRxns.rxnNames(idx) = [];
0165             newRxns.equations(idx) = [];
0166             newRxns.grRules(idx) = [];
0167         end
0168     end
0169 
0170     newRxns.lb = zeros(length(newRxns.rxns),1);
0171     newRxns.ub = zeros(length(newRxns.rxns),1) + 1000;
0172 
0173     model = addRxns(model, newRxns, 3, '', true, true);
0174 
0175     % 4. Add protein usage reactions.
0176     usageRxns.rxns            = strcat('usage_',proteinMets.mets);
0177     usageRxns.rxnNames        = usageRxns.rxns;
0178     usageRxns.mets            = cell(numel(usageRxns.rxns),1);
0179     usageRxns.stoichCoeffs    = cell(numel(usageRxns.rxns),1);
0180     for i=1:numel(usageRxns.mets)
0181         usageRxns.mets{i}         = {proteinMets.mets{i}, 'prot_pool'};
0182         usageRxns.stoichCoeffs{i} = [-1,1];
0183     end
0184     usageRxns.lb              = zeros(numel(usageRxns.rxns),1) - 1000;
0185     usageRxns.ub              = zeros(numel(usageRxns.rxns),1);
0186     usageRxns.rev             = ones(numel(usageRxns.rxns),1);
0187     usageRxns.grRules         = newEnzymes.genes;
0188     model = addRxns(model,usageRxns);
0189 
0190     % 5. Update .ec field with enzyme data
0191     model.ec.enzymes    = [model.ec.enzymes; newEnzymes.enzymes(:)];
0192     model.ec.genes      = [model.ec.genes; newEnzymes.genes(:)];
0193     model.ec.mw         = [model.ec.mw; newEnzymes.mw(:)];
0194     model.ec.sequence   = [model.ec.sequence; repmat({''},numel(newEnzymes.enzymes),1)]; 
0195     if isfield(model.ec,'concs')
0196         model.ec.concs  = [model.ec.concs; NaN(numel(newEnzymes.enzymes),1)];
0197     end
0198 
0199     % 6. Expand the enzyme rxns matrix
0200     rxnWithRule = find(~cellfun('isempty',newRxns.grRules));
0201 
0202     model.ec.rxnEnzMat =  [model.ec.rxnEnzMat, zeros(length(model.ec.rxns), length(newEnzymes.enzymes))]; % add n new enzymes
0203     model.ec.rxnEnzMat =  [model.ec.rxnEnzMat; zeros(length(rxnWithRule), length(model.ec.enzymes))]; % add n new rxns
0204 
0205     % 7. Update .ec field with rxn data
0206 
0207     % Get reactions to add
0208     numRxns = length(model.ec.rxns);
0209     for i = 1:numel(rxnWithRule)
0210         rxnIdx = rxnWithRule(i);
0211 
0212         model.ec.rxns(end+1)    = newRxns.rxns(rxnIdx);
0213         model.ec.kcat(end+1)    = 0;
0214         model.ec.source(end+1)  = {''};
0215         model.ec.notes(end+1)   = {''};
0216         model.ec.eccodes(end+1) = {''};
0217 
0218         % Get rules
0219         genes = newRxns.grRules(rxnIdx);
0220         genes = strtrim(erase(split(genes,'and'), {'(' ')'}));
0221         [~,enzIdx] = ismember(genes, model.ec.genes);
0222 
0223         % Update the enzyme rxns matrix
0224         model.ec.rxnEnzMat(numRxns+i, enzIdx) = 1;
0225     end
0226     
0227     rxnsAdded = newRxns.rxns(:);
0228 end
0229 end