0001 function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,...
0002 deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, reservedRxns, suppressWarnings)
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037 if nargin<2
0038 deleteUnconstrained=true;
0039 end
0040 if nargin<3
0041 deleteDuplicates=false;
0042 end
0043 if nargin<4
0044 deleteZeroInterval=false;
0045 end
0046 if nargin<5
0047 deleteInaccessible=false;
0048 end
0049 if nargin<6
0050 deleteMinMax=false;
0051 end
0052 if nargin<7
0053 groupLinear=false;
0054 end
0055 if nargin<8
0056 reservedRxns=[];
0057 end
0058 if nargin<9
0059 suppressWarnings=false;
0060 end
0061
0062 reducedModel=model;
0063 deletedReactions={};
0064 deletedMetabolites={};
0065
0066 if deleteUnconstrained==true
0067 if isfield(reducedModel,'unconstrained')
0068
0069 deletedMetabolites=reducedModel.mets(reducedModel.unconstrained~=0);
0070 reducedModel=removeMets(reducedModel,reducedModel.unconstrained~=0);
0071 reducedModel=rmfield(reducedModel,'unconstrained');
0072 end
0073 end
0074
0075 if deleteDuplicates==true
0076
0077
0078
0079 [reducedModel rxnsToDelete]=contractModel(reducedModel);
0080 deletedReactions=[deletedReactions; rxnsToDelete];
0081 end
0082
0083 if deleteZeroInterval==true
0084
0085 zeroIntervalReactions=and(reducedModel.lb==0,reducedModel.ub==0);
0086
0087 rxnsToDelete=setdiff(reducedModel.rxns(zeroIntervalReactions),reservedRxns);
0088 deletedReactions=[deletedReactions; rxnsToDelete];
0089
0090
0091 reducedModel=removeRxns(reducedModel,rxnsToDelete);
0092
0093
0094 notInUse=sum(reducedModel.S~=0,2)==0;
0095 deletedMetabolites=[deletedMetabolites;reducedModel.mets(notInUse)];
0096
0097
0098 reducedModel=removeMets(reducedModel,notInUse);
0099 end
0100
0101 if deleteInaccessible==true
0102
0103
0104
0105 if isfield(reducedModel,'unconstrained') && suppressWarnings==false
0106 dispEM('Removing dead-end reactions before removing exchange metabolites',false);
0107 end
0108
0109 while true
0110
0111
0112
0113
0114
0115
0116
0117
0118 in=any(reducedModel.b<0,2);
0119 out=any(reducedModel.b>0,2);
0120 I=speye(numel(reducedModel.mets));
0121 revS=[reducedModel.S,reducedModel.S(:,reducedModel.rev~=0)*-1 I(:,in) I(:,out)*-1];
0122
0123 metUsage=sum(abs(revS')>0);
0124 onlyProducts=sum(revS'>0) == metUsage;
0125 onlyReactants=sum(revS'<0) == metUsage;
0126
0127
0128
0129
0130 notInUse=onlyProducts | onlyReactants | (sum(abs(reducedModel.S')>0)<=1 & (~in & ~out)');
0131 deletedRxn=false;
0132 if any(notInUse)
0133
0134 rxnsNotInUse=sum(abs(reducedModel.S(notInUse,:))>0,1)>0;
0135 rxnsToDelete=setdiff(reducedModel.rxns(rxnsNotInUse),reservedRxns);
0136 deletedReactions=[deletedReactions; rxnsToDelete];
0137
0138
0139 reducedModel=removeRxns(reducedModel,rxnsToDelete);
0140
0141
0142
0143 notInUse=sum(reducedModel.S~=0,2)==0;
0144 deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0145 reducedModel=removeMets(reducedModel,notInUse);
0146
0147
0148
0149 if ~isempty(rxnsToDelete)
0150 deletedRxn=true;
0151 end
0152 else
0153 break;
0154 end
0155 if deletedRxn==false
0156 break;
0157 end
0158 end
0159 end
0160
0161 if deleteMinMax==true
0162
0163
0164 I=~haveFlux(reducedModel);
0165
0166
0167 rxnsToDelete=setdiff(reducedModel.rxns(I),reservedRxns);
0168 deletedReactions=[deletedReactions; rxnsToDelete];
0169 reducedModel=removeRxns(reducedModel,rxnsToDelete);
0170
0171
0172 notInUse=sum(reducedModel.S~=0,2)==0;
0173 deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0174 reducedModel=removeMets(reducedModel,notInUse);
0175 end
0176
0177
0178 if groupLinear==true
0179 if ~any(reducedModel.rev)
0180 if suppressWarnings==false
0181 fprintf('NOTE: You have chosen to group linear reactions. This option does not keep track of gene/reaction associations when reactions are merged. Deleting all gene information\n');
0182 end
0183 bannedIndexes=getIndexes(reducedModel,reservedRxns,'rxns');
0184 while 1
0185
0186
0187 singleNegative=find(sum(reducedModel.S'<0)==1);
0188 singlePositive=find(sum(reducedModel.S'>0)==1);
0189
0190
0191 common=intersect(singleNegative,singlePositive);
0192
0193
0194 mergedSomeRxn=false;
0195 for i=1:numel(common)
0196 involvedRxns=find(reducedModel.S(common(i),:));
0197
0198
0199
0200
0201
0202 if isempty(intersect(bannedIndexes,involvedRxns)) && any(involvedRxns)
0203
0204
0205
0206 stoichRatio=abs(reducedModel.S(common(i),involvedRxns(1))/reducedModel.S(common(i),involvedRxns(2)));
0207 reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio;
0208 reducedModel.S(common(i),involvedRxns(1))=0;
0209
0210
0211 reducedModel.rxns{involvedRxns(1)}=[reducedModel.rxns{involvedRxns(1)} '_' reducedModel.rxns{involvedRxns(2)}];
0212 reducedModel.rxnNames{involvedRxns(1)}=reducedModel.rxns{involvedRxns(1)};
0213
0214
0215
0216 reducedModel.S(:,involvedRxns(2))=0;
0217 mergedSomeRxn=true;
0218 end
0219 end
0220
0221
0222 if mergedSomeRxn==false
0223 break;
0224 end
0225 end
0226
0227
0228 I=sum(reducedModel.S~=0)==0;
0229
0230
0231 deletedReactions=[deletedReactions; reducedModel.rxns(I)];
0232 reducedModel=removeRxns(reducedModel,find(I));
0233
0234
0235 notInUse=sum(reducedModel.S~=0,2)==0;
0236 deletedMetabolites=[deletedMetabolites; reducedModel.mets(notInUse)];
0237 reducedModel=removeMets(reducedModel,notInUse);
0238
0239 reducedModel.genes={};
0240 reducedModel.rxnGeneMat=sparse(numel(reducedModel.rxns),0);
0241 reducedModel.grRules(:)={''};
0242
0243 if isfield(reducedModel,'geneShortNames')
0244 reducedModel.geneShortNames={};
0245 end
0246 if isfield(reducedModel,'geneMiriams')
0247 reducedModel.geneMiriams={};
0248 end
0249 if isfield(reducedModel,'geneComps')
0250 reducedModel.geneComps=[];
0251 end
0252 else
0253 dispEM('Cannot group reactions for a model that has reversible reactions');
0254 end
0255 end
0256
0257 end