removeGenes

 removeGenes
   Deletes a set of genes from a model

   model                   a model structure
   genesToRemove           either a cell array of gene IDs, a logical vector
                           with the same number of elements as genes in the model,
                           or a vector of indexes to remove
   removeUnusedMets        remove metabolites that are no longer in use (opt, default
                           false)
   removeBlockedRxns       remove reactions that get blocked after deleting the genes
                           (opt, default false)
   standardizeRules        format gene rules to be compliant with standard format
                           (opt, default true)

   reducedModel            an updated model structure

   Usage: reducedModel = removeGenes(model,genesToRemove,removeUnusedMets,removeBlockedRxns)

0001 function reducedModel = removeGenes(model,genesToRemove,removeUnusedMets,removeBlockedRxns,standardizeRules)
0002 % removeGenes
0003 %   Deletes a set of genes from a model
0004 %
0005 %   model                   a model structure
0006 %   genesToRemove           either a cell array of gene IDs, a logical vector
0007 %                           with the same number of elements as genes in the model,
0008 %                           or a vector of indexes to remove
0009 %   removeUnusedMets        remove metabolites that are no longer in use (opt, default
0010 %                           false)
0011 %   removeBlockedRxns       remove reactions that get blocked after deleting the genes
0012 %                           (opt, default false)
0013 %   standardizeRules        format gene rules to be compliant with standard format
0014 %                           (opt, default true)
0015 %
0016 %   reducedModel            an updated model structure
0017 %
0018 %   Usage: reducedModel = removeGenes(model,genesToRemove,removeUnusedMets,removeBlockedRxns)
0019 
0020 if nargin<3
0021     removeUnusedMets = false;
0022 end
0023 if nargin<4
0024     removeBlockedRxns = false;
0025 end
0026 if nargin<5
0027     standardizeRules = true;
0028 end
0029 %Format grRules and rxnGeneMatrix:
0030 if standardizeRules
0031     [grRules,rxnGeneMat,toCheck] = standardizeGrRules(model,true);
0032     model.grRules    = grRules;
0033     model.rxnGeneMat = rxnGeneMat;
0034 else
0035     toCheck    = [];
0036     rxnGeneMat = model.rxnGeneMat;
0037 end
0038 reducedModel = model;
0039 %Only remove genes that are actually in the model
0040 if ~(islogical(genesToRemove) || isnumeric(genesToRemove))
0041     genesToRemove=convertCharArray(genesToRemove);
0042     genesToRemove=genesToRemove(ismember(genesToRemove,model.genes));
0043 end
0044 if ~isempty(genesToRemove)
0045     indexesToRemove = getIndexes(model,genesToRemove,'genes');
0046     if ~isempty(indexesToRemove)
0047         %Make 0 corresponding columns from rxnGeneMat:
0048         reducedModel.rxnGeneMat(:,indexesToRemove) = 0;        
0049         genes = model.genes(indexesToRemove);
0050         %Check if conflicting grRules contain any of the genes to remove
0051         if ~isempty(toCheck)
0052             for i=1:numel(toCheck)
0053                 index   = toCheck(i);
0054                 gIdxs   = find(rxnGeneMat(index,:));
0055                 g2check = model.genes(gIdxs);
0056                 if any(ismember(g2check,genes))
0057                     warning(['Conflicting grRule #' num2str(index) ' (' model.rxns{index} ') contains at least one of the genes to be removed, this grRule will be bypassed in order to avoid logical errors'])
0058                 end
0059             end
0060         end
0061         canCarryFlux = true(size(model.rxns));       
0062         %Loop through genes and adapt rxns:
0063         for i = 1:length(genes)
0064             %Find all reactions for this gene and loop through them:
0065             geneRxns = find(rxnGeneMat(:,indexesToRemove(i)));
0066             if ~isempty(geneRxns)
0067                 for j = 1:numel(geneRxns)
0068                     index  = geneRxns(j);
0069                     grRule = reducedModel.grRules{index};
0070                     ruleGenes = reducedModel.genes(logical(rxnGeneMat(index,:)));
0071                     if ~ismember(index,toCheck) && canCarryFlux(index) && ~isempty(grRule)
0072                         %Check if rxn can carry flux without this gene:
0073                         canCarryFlux(index) = canRxnCarryFlux(ruleGenes,grRule,genes{i});
0074                         %Adapt gene rule & gene matrix:
0075                         grRule = removeGeneFromRule(grRule,genes{i});
0076                         reducedModel.grRules{index} = grRule;
0077                     end
0078                 end
0079             end
0080         end
0081         %Delete or block the reactions that cannot carry flux:
0082         if removeBlockedRxns
0083             rxnsToRemove = reducedModel.rxns(~canCarryFlux);
0084             reducedModel = removeReactions(reducedModel,rxnsToRemove,removeUnusedMets,true);
0085         else
0086             reducedModel = removeReactions(reducedModel,[],removeUnusedMets,true);
0087             reducedModel.lb(~canCarryFlux) = 0;
0088             reducedModel.ub(~canCarryFlux) = 0;
0089         end
0090     end
0091 end
0092 %Format grRules and rxnGeneMatrix after all modifications
0093 if standardizeRules
0094     [grRules,rxnGeneMat] = standardizeGrRules(reducedModel,true);
0095     reducedModel.grRules = grRules;
0096     reducedModel.rxnGeneMat = rxnGeneMat;
0097 end
0098 end
0099 
0100 function canIt = canRxnCarryFlux(ruleGenes,geneRule,geneToRemove)
0101 %This function converts a gene rule to a logical statement, and then
0102 %asseses if the rule is true (i.e. rxn can still carry flux) or not (cannot
0103 %carry flux).
0104 geneRule = [' ', geneRule, ' '];
0105 for i = 1:length(ruleGenes)
0106     if strcmp(ruleGenes{i},geneToRemove)
0107         geneRule = strrep(geneRule,[' ' ruleGenes{i} ' '],' false ');
0108         geneRule = strrep(geneRule,['(' ruleGenes{i} ' '],'(false ');
0109         geneRule = strrep(geneRule,[' ' ruleGenes{i} ')'],' false)');
0110     else
0111         geneRule = strrep(geneRule,[' ' ruleGenes{i} ' '],' true ');
0112         geneRule = strrep(geneRule,['(' ruleGenes{i} ' '],'(true ');
0113         geneRule = strrep(geneRule,[' ' ruleGenes{i} ')'],' true)');
0114     end
0115 end
0116 geneRule = strtrim(geneRule);
0117 geneRule = strrep(geneRule,'and','&&');
0118 geneRule = strrep(geneRule,'or','||');
0119 canIt    = eval(geneRule);
0120 end
0121 
0122 function geneRule = removeGeneFromRule(geneRule,geneToRemove)
0123 %This function receives a standard gene rule and it returns it without the
0124 %chosen gene.
0125 geneSets = strsplit(geneRule,' or ');
0126 hasGene  = ~cellfun(@isempty,strfind(geneSets,geneToRemove));
0127 geneSets = geneSets(~hasGene);
0128 geneRule = strjoin(geneSets,' or ');
0129 end