removeLowScoreGenes Remove low-scoring genes from model.
This function removes genes from a model based on their scores, a step
used by the tINIT package. The function recognizes and differentiates
between isozymes and subunits of an enzyme complex. Genes are removed
from each grRule, subject to the following conditions:
1) At least one gene must remain associated with the reaction
2) Genes involved in a complex (joined by ANDs) are not removed
Usage:
[newModel,remGenes] = removeLowScoreGenes(model,geneScores,complexScoring);
Inputs:
model Model structure from which genes are to be removed.
geneScores A vector of scores associated with the model genes.
Genes with a positive score will remain in the model,
whereas genes with a negative score will try to be
removed.
If all genes associated with a reaction have a negative
score, then the least-negative gene will remain; if
there is a tie, one will be selected at random.
If a negative-scoring gene is a subunit in a complex,
it will not be removed; however, the entire complex may
be removed. See the following example cases:
Original: G1 or (G2 and G3 and G4)
Negative: G1, G2
New: G2 and G3 and G4
Original: G1 or (G2 and G3) or (G4 and G5)
Negative: G1, G2
New: G4 and G5 [using default complexScoring]
Original: (G1 and (G2 or G3) and G4)
Negative: G2
New: (G1 and G3 and G4)
isozymeScoring Method used to combine the scores of multiple isozymes;
'min', 'max', 'median', or 'average'.
(optional, default 'max')
complexScoring Method used to combine the scores of enzyme complex
subunits: 'min', 'max', 'median', or 'average'.
(optional, default 'min')
Outputs:
newModel Model with updated gene, grRules, and rxnGeneMat fields
after attemping to remove negative-score genes.
remGenes A list of negative-score genes that were fully removed
from the model. Negative-score genes that were removed
from some grRules but not all will not be included in
this list.0001 function [newModel,remGenes] = removeLowScoreGenes(model,geneScores,isozymeScoring,complexScoring) 0002 %removeLowScoreGenes Remove low-scoring genes from model. 0003 % 0004 % This function removes genes from a model based on their scores, a step 0005 % used by the tINIT package. The function recognizes and differentiates 0006 % between isozymes and subunits of an enzyme complex. Genes are removed 0007 % from each grRule, subject to the following conditions: 0008 % 1) At least one gene must remain associated with the reaction 0009 % 2) Genes involved in a complex (joined by ANDs) are not removed 0010 % 0011 % Usage: 0012 % 0013 % [newModel,remGenes] = removeLowScoreGenes(model,geneScores,complexScoring); 0014 % 0015 % Inputs: 0016 % 0017 % model Model structure from which genes are to be removed. 0018 % 0019 % geneScores A vector of scores associated with the model genes. 0020 % Genes with a positive score will remain in the model, 0021 % whereas genes with a negative score will try to be 0022 % removed. 0023 % 0024 % If all genes associated with a reaction have a negative 0025 % score, then the least-negative gene will remain; if 0026 % there is a tie, one will be selected at random. 0027 % 0028 % If a negative-scoring gene is a subunit in a complex, 0029 % it will not be removed; however, the entire complex may 0030 % be removed. See the following example cases: 0031 % 0032 % Original: G1 or (G2 and G3 and G4) 0033 % Negative: G1, G2 0034 % New: G2 and G3 and G4 0035 % 0036 % Original: G1 or (G2 and G3) or (G4 and G5) 0037 % Negative: G1, G2 0038 % New: G4 and G5 [using default complexScoring] 0039 % 0040 % Original: (G1 and (G2 or G3) and G4) 0041 % Negative: G2 0042 % New: (G1 and G3 and G4) 0043 % 0044 % isozymeScoring Method used to combine the scores of multiple isozymes; 0045 % 'min', 'max', 'median', or 'average'. 0046 % (optional, default 'max') 0047 % 0048 % complexScoring Method used to combine the scores of enzyme complex 0049 % subunits: 'min', 'max', 'median', or 'average'. 0050 % (optional, default 'min') 0051 % 0052 % Outputs: 0053 % 0054 % newModel Model with updated gene, grRules, and rxnGeneMat fields 0055 % after attemping to remove negative-score genes. 0056 % 0057 % remGenes A list of negative-score genes that were fully removed 0058 % from the model. Negative-score genes that were removed 0059 % from some grRules but not all will not be included in 0060 % this list. 0061 % 0062 0063 0064 if nargin < 3 || isempty(isozymeScoring) 0065 isozymeScoring = 'max'; 0066 end 0067 if nargin < 4 0068 complexScoring = 'min'; 0069 end 0070 0071 if ~isequal(size(model.genes),size(geneScores)) 0072 error('The dimensions of model genes and geneScores do not match.'); 0073 end 0074 0075 % convert logical operators to symbols 0076 grRules = model.grRules; 0077 if any(contains(grRules,{'&','|'})) 0078 symbolicRules = true; 0079 else 0080 symbolicRules = false; 0081 end 0082 grRules = regexprep(grRules,' and ',' & '); 0083 grRules = regexprep(grRules,' or ',' | '); 0084 0085 % get unique list of grRules 0086 [uRules,~,rule_ind] = unique(grRules); 0087 0088 % iterate through each rule 0089 newRules = uRules; %initialize newRules 0090 for i = 1:numel(uRules) 0091 if isempty(uRules{i}) || ~contains(uRules{i},'|') 0092 continue 0093 elseif contains(uRules{i},'&') 0094 newRules{i} = processComplexRule(uRules{i},model.genes,geneScores,isozymeScoring,complexScoring); 0095 else 0096 newRules{i} = processSimpleRule(uRules{i},model.genes,geneScores,isozymeScoring,complexScoring); 0097 end 0098 end 0099 0100 % re-map unique rules to model 0101 newModel = model; 0102 newModel.grRules = newRules(rule_ind); 0103 0104 % restore original logical operator formatting if it was changed 0105 if ~symbolicRules 0106 newModel.grRules = regexprep(newModel.grRules, ' \| ', ' or '); 0107 newModel.grRules = regexprep(newModel.grRules, ' & ', ' and '); 0108 end 0109 0110 % regenerate "genes" and "rxnGeneMat" model fields 0111 [genes,rxnGeneMat] = getGenesFromGrRules(newModel.grRules); 0112 0113 % determine which of the original genes were removed 0114 remInd = ~ismember(model.genes,genes); 0115 remGenes = model.genes(remInd); 0116 0117 % Keep the retained genes in their original order rather than the sorted 0118 % order returned by getGenesFromGrRules. Gene removal never introduces new 0119 % genes, so the remaining genes are model.genes(~remInd). Preserving this 0120 % order ensures the gene-associated fields trimmed below (geneShortNames, 0121 % proteins, etc.), which are indexed by remInd, stay aligned with genes. 0122 newModel.genes = model.genes(~remInd); 0123 [~,reorderInd] = ismember(newModel.genes,genes); 0124 newModel.rxnGeneMat = rxnGeneMat(:,reorderInd); 0125 0126 if isfield(newModel,'geneShortNames') 0127 newModel.geneShortNames(remInd) = []; 0128 end 0129 if isfield(newModel,'proteins') 0130 newModel.proteins(remInd) = []; 0131 end 0132 if isfield(newModel,'geneMiriams') 0133 newModel.geneMiriams(remInd) = []; 0134 end 0135 if isfield(newModel,'geneFrom') 0136 newModel.geneFrom(remInd) = []; 0137 end 0138 if isfield(newModel,'geneComps') 0139 newModel.geneComps(remInd) = []; 0140 end 0141 0142 0143 end 0144 0145 0146 0147 function [updatedRule,rScore] = processSimpleRule(rule,genes,gScores,isozymeScoring,complexScoring) 0148 % Either score or modify a reaction gene rule containig only ANDs or ORs. 0149 % 0150 % If the rule contains an enzyme complex (all ANDs), the complex will be 0151 % scored based on the score of its subunits. Subunits without a score (NaN) 0152 % will be excluded from the score calculation. 0153 % 0154 % If the rule contains only isozymes (all ORs), the negative-score genes 0155 % will be removed from the rule. Isozymes without a score (NaN) will not be 0156 % removed from the rule. The resuling rule will then be scored. 0157 0158 0159 % get IDs and indices of genes involved in rule 0160 ruleGenes = unique(regexp(rule,'[^&|\(\) ]+','match')); 0161 [~,geneInd] = ismember(ruleGenes,genes); 0162 0163 % rules with one or no genes remain unchanged 0164 if numel(ruleGenes) < 2 0165 rScore = gScores(geneInd); 0166 updatedRule = rule; 0167 return 0168 end 0169 0170 if ~contains(rule,'&') % rule contains isozymes 0171 0172 scoreMethod = isozymeScoring; 0173 negInd = gScores(geneInd) < 0; % NaNs will return false here 0174 if all(negInd) 0175 % get the least negative gene, adding a small random value to avoid a tie 0176 [~,maxInd] = max(gScores(geneInd) + rand(size(geneInd))*(1e-8)); 0177 updatedRule = ruleGenes{maxInd}; 0178 elseif sum(~negInd) == 1 0179 updatedRule = ruleGenes{~negInd}; 0180 else 0181 updatedRule = strjoin(ruleGenes(~negInd),' | '); 0182 if startsWith(rule,'(') 0183 updatedRule = ['(',updatedRule,')']; 0184 end 0185 end 0186 0187 % update ruleGenes and their indices 0188 ruleGenes = unique(regexp(updatedRule,'[^&|\(\) ]+','match')); 0189 [~,geneInd] = ismember(ruleGenes,genes); 0190 0191 elseif ~contains(rule,'|') % rule contains enzyme complex 0192 scoreMethod = complexScoring; 0193 updatedRule = rule; 0194 else 0195 error('This function cannot handle rules with both "OR" and "AND" expressions.'); 0196 end 0197 0198 % score rule 0199 switch lower(scoreMethod) 0200 case 'min' 0201 rScore = min(gScores(geneInd),[],'omitnan'); 0202 case 'max' 0203 rScore = max(gScores(geneInd),[],'omitnan'); 0204 case 'median' 0205 rScore = median(gScores(geneInd),'omitnan'); 0206 case 'average' 0207 rScore = mean(gScores(geneInd),'omitnan'); 0208 end 0209 0210 end 0211 0212 0213 0214 function updatedRule = processComplexRule(rule,genes,gScores,isozymeScoring,complexScoring) 0215 % Update reactions containing both AND and OR expressions. 0216 % 0217 % Negative-score genes will be removed if they are isozymic, whereas they 0218 % will not be removed if they are part of an enzyme complex. However, if 0219 % the enzyme complex has a negative score, the entire complex will be 0220 % removed, as long as it is not the only remaining element in the rule. 0221 0222 0223 % Specify phrases to search for in the grRule. These phrases will find 0224 % genes grouped by all ANDs (first phrase) or all ORs (second phrase). 0225 search_phrases = {'\([^&|\(\) ]+( & [^&|\(\) ]+)+\)', '\([^&|\(\) ]+( \| [^&|\(\) ]+)+\)'}; 0226 0227 % initialize some variables 0228 subsets = {}; % subsets are groups of genes grouped by all ANDs or all ORs 0229 c = 1; % counter to keep track of the group (subset) number 0230 r_orig = rule; % record original rule to determine when it stops changing 0231 for k = 1:100 % iterate some arbitrarily high number of times 0232 for j = 1:length(search_phrases) 0233 new_subset = regexp(rule,search_phrases{j},'match')'; % extract subsets 0234 if ~isempty(new_subset) 0235 subsets = [subsets; new_subset]; % append to list of subsets 0236 subset_nums = arrayfun(@num2str,(c:length(subsets))','UniformOutput',false); % get group numbers to be assigned to the new subsets, and convert to strings 0237 rule = regexprep(rule,search_phrases{j},strcat('#',subset_nums,'#'),'once'); % replace the subsets in the expression with their group numbers (enclosed by "#"s) 0238 c = c + length(new_subset); 0239 end 0240 end 0241 if isequal(rule,r_orig) 0242 break; % stop iterating when rule stops changing 0243 else 0244 r_orig = rule; 0245 end 0246 end 0247 subsets{end+1} = rule; % add final state of rule as the last subset 0248 0249 % score and update each subset, and append to gene list and gene scores 0250 for i = 1:numel(subsets) 0251 [subsets{i},subset_score] = processSimpleRule(subsets{i},genes,gScores,isozymeScoring,complexScoring); 0252 gScores = [gScores; subset_score]; 0253 genes = [genes; {strcat('#',num2str(i),'#')}]; 0254 end 0255 0256 % reconstruct the rule from its updated subsets 0257 updatedRule = subsets{end}; 0258 for i = c-1:-1:1 0259 updatedRule = regexprep(updatedRule,strcat('#',num2str(i),'#'),subsets{i}); 0260 end 0261 0262 end 0263 0264 0265