mergeCompartments
Merge all compartments in a model
model a model structure
keepUnconstrained keep metabolites that are unconstrained in a
'unconstrained' compartment. If these are merged the
exchange reactions will most often be deleted (optional,
default false)
deleteRxnsWithOneMet delete reactions with only one metabolite. These
reactions come from reactions such as A[c] + B[c]
=> A[m]. In some models hydrogen is balanced around
each membrane with reactions like this (optional,
default false)
distReverse distinguish reactions with same metabolites but
different reversibility as different reactions
(optional, default true)
model a model with all reactions located to one compartment
deletedRxns reactions that were deleted because of only
having one metabolite after merging
duplicateRxns identical reactions that occurred in different
compartments and were deleted because they turned
to be duplicated after merging
Merges all compartments into one 's' compartment (for 'System'). This can
be useful for example to ensure that there are metabolic capabilities to
synthesize all metabolites.
NOTE: If the metabolite IDs reflect the compartment that they are in
the IDs may no longer be representative.
Usage: [model, deletedRxns, duplicateRxns]=mergeCompartments(model,keepUnconstrained,deleteRxnsWithOneMet,distReverse)0001 function [model, deletedRxns, duplicateRxns]=mergeCompartments(model,keepUnconstrained,deleteRxnsWithOneMet,distReverse) 0002 % mergeCompartments 0003 % Merge all compartments in a model 0004 % 0005 % model a model structure 0006 % keepUnconstrained keep metabolites that are unconstrained in a 0007 % 'unconstrained' compartment. If these are merged the 0008 % exchange reactions will most often be deleted (optional, 0009 % default false) 0010 % deleteRxnsWithOneMet delete reactions with only one metabolite. These 0011 % reactions come from reactions such as A[c] + B[c] 0012 % => A[m]. In some models hydrogen is balanced around 0013 % each membrane with reactions like this (optional, 0014 % default false) 0015 % distReverse distinguish reactions with same metabolites but 0016 % different reversibility as different reactions 0017 % (optional, default true) 0018 % 0019 % model a model with all reactions located to one compartment 0020 % deletedRxns reactions that were deleted because of only 0021 % having one metabolite after merging 0022 % duplicateRxns identical reactions that occurred in different 0023 % compartments and were deleted because they turned 0024 % to be duplicated after merging 0025 % 0026 % Merges all compartments into one 's' compartment (for 'System'). This can 0027 % be useful for example to ensure that there are metabolic capabilities to 0028 % synthesize all metabolites. 0029 % 0030 % NOTE: If the metabolite IDs reflect the compartment that they are in 0031 % the IDs may no longer be representative. 0032 % 0033 % Usage: [model, deletedRxns, duplicateRxns]=mergeCompartments(model,keepUnconstrained,deleteRxnsWithOneMet,distReverse) 0034 0035 if nargin<2 0036 keepUnconstrained=false; 0037 end 0038 if nargin<3 0039 deleteRxnsWithOneMet=false; 0040 end 0041 if nargin<4 0042 distReverse=true; 0043 end 0044 0045 if ~isfield(model,'unconstrained') 0046 keepUnconstrained=false; 0047 end 0048 0049 %Keep track of which reactions only contained one metabolite before the 0050 %merging as they are probable exchange reactions. 0051 if deleteRxnsWithOneMet==true 0052 reservedRxns=model.rxns(sum(model.S~=0)==1); 0053 if ~isempty(reservedRxns) && isfield(model,'unconstrained') 0054 %If there is no unconstrained field these reactions are probably 0055 %exchange reactions and shall be kept. If not then print a warning 0056 EM='There are reactions with only one metabolite. Cannot determine whether they are exchange reactions since there is no unconstrained field'; 0057 dispEM(EM,false); 0058 end 0059 end 0060 0061 %Loop through each metabolite, and if it is not unconstrained then change 0062 %the S matrix to use the metabolite with the lowest index in model.comps 0063 %instead 0064 uNames=unique(model.metNames); 0065 for i=1:numel(uNames) 0066 %Find all metabolites with this name.. 0067 I=ismember(model.metNames,uNames(i)); 0068 0069 %Find the first of those that is not unconstrained. This is the one 0070 %that the other "un-unconstrained" should be changed to. 0071 if keepUnconstrained==true 0072 mergeTo=find(I & model.unconstrained==false,1); 0073 0074 %This could happen if there is only one metabolite and it is 0075 %unconstrained 0076 if isempty(mergeTo) 0077 continue; 0078 end 0079 else 0080 mergeTo=find(I,1); 0081 end 0082 I(mergeTo)=false; %Do not do anything for itself 0083 I=find(I); 0084 0085 %Go through each of the metabolites with this name and update them to 0086 %mergeTo 0087 for j=1:numel(I) 0088 if keepUnconstrained==true && model.unconstrained(I(j))==true 0089 continue; 0090 end 0091 %Update the S matrix 0092 model.S(mergeTo,:)=model.S(mergeTo,:)+model.S(I(j),:); 0093 model.S(I(j),:)=0; 0094 end 0095 end 0096 0097 %Remove all metabolites that are no longer used (with a bit of a trick) 0098 model=removeReactions(model,{},true); 0099 0100 %Update all mets so that they are in compartment "s" with id "1" 0101 model.compNames={'System'}; 0102 model.comps={'s'}; 0103 model.compOutside={''}; 0104 model.metComps=ones(numel(model.mets),1); 0105 0106 %Add exchange mets to another compartment "b" with id "2" 0107 if keepUnconstrained==true 0108 model.compNames=[model.compNames; {'Unconstrained'}]; 0109 model.comps=[model.comps; {'b'}]; 0110 model.compOutside=[model.compOutside; {'s'}]; 0111 model.metComps(model.unconstrained~=0)=2; 0112 end 0113 0114 %Transport reactions on the form A <=> B will have been deleted by the 0115 %merging. Remove those reactions 0116 model=removeMets(model,{},false,true,true); 0117 0118 if deleteRxnsWithOneMet==true 0119 I=model.rxns(sum(model.S~=0)==1); 0120 0121 %Delete the reactions that contain only one metabolite after the 0122 %merging but not before 0123 deletedRxns=setdiff(I,reservedRxns); 0124 model=removeReactions(model,deletedRxns,true,true); 0125 else 0126 deletedRxns={}; 0127 end 0128 0129 %And then finally merge the identical reactions 0130 [model, ~, duplicateRxns]=contractModel(model,distReverse); 0131 end