diff --git a/INIT/ftINIT.m b/INIT/ftINIT.m index ab5f3e030..9dba76ad1 100644 --- a/INIT/ftINIT.m +++ b/INIT/ftINIT.m @@ -135,7 +135,7 @@ %Get the matching mets metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers - (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers + (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' ~= 0); %reversible reactions that are consumers %convert the production rxns from refModel to minModel prepData.groupIds; [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); diff --git a/INIT/mergeLinear.m b/INIT/mergeLinear.m index 1378353b2..114d7228d 100644 --- a/INIT/mergeLinear.m +++ b/INIT/mergeLinear.m @@ -62,8 +62,12 @@ for i=1:numel(twoNonZero) involvedRxns=find(reducedModel.S(twoNonZero(i),:)); %Check that we can have one positive and one negative - pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns)); - neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns)); + %rev is compared rather than used as-is: the other operand is a + %sparse logical, and MATLAB refuses to combine one with an integer + %array ("Sparse integer array arithmetic operations are not + %supported"). A model built outside RAVEN can carry an integer rev. + pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0); + neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0); %Check so that one or both of the reactions haven't been merged diff --git a/core/replaceMets.m b/core/replaceMets.m index 9301d8628..9cfa0804c 100755 --- a/core/replaceMets.m +++ b/core/replaceMets.m @@ -20,7 +20,7 @@ % Output: % model model structure with selected metabolites replaced % removedRxns identifiers of duplicate reactions that were removed -% idxDuplRxns index of removedRxns in original model +% idxDuplRxns indices of removedRxns in the input model % % Note: This function is useful when the model contains both 'oxygen' and % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, @@ -59,10 +59,11 @@ error('The to-be-replaced metabolite cannot be found in the model.'); end -rxnsWithMet = find(model.S(metIdx,:)); +[~, rxnsWithMet] = find(model.S(metIdx,:)); +rxnsWithMet = unique(rxnsWithMet); if verbose==true - fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') - fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) + fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n') + fprintf('%s\n',model.rxns{rxnsWithMet}) fprintf('\n') end @@ -88,8 +89,11 @@ idxDelete=[]; if identifiers + if any(ismember(metIdx,repIdx)) + error('The to-be-replaced metabolite and the replacement metabolite are the same.'); + end originalStoch = model.S(metIdx,rxnsWithMet); - model.S(repIdx,rxnsWithMet) = originalStoch; + model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch; model.S(metIdx,rxnsWithMet) = 0; idxDelete = metIdx; else @@ -99,54 +103,38 @@ % the to-be-replace metabolite ID deleted. % Build list of metaboliteName[compartment] - metCompsN =cellstr(num2str(model.metComps)); - map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); - metCompsN = map.values(metCompsN); - metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); + comps = model.comps(:); + metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']'); - for i = 1:length(repIdx) - metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); - if length(metCompsNidx)>1 - for j = 2:length(metCompsNidx) - model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); - idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete - end + % Within each metaboliteName[compartment] group, the first replacement + % metabolite is kept and the other members are merged into it. + groups = unique(metCompsN(repIdx),'stable'); + for i = 1:numel(groups) + metCompsNidx = find(strcmp(groups{i}, metCompsN)); + keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1)); + others = metCompsNidx(metCompsNidx ~= keepIdx); + if ~isempty(others) + model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1); + idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete end end end - if ~isempty(idxDelete) - model.S(idxDelete,:) =[]; - model.mets(idxDelete) = []; - model.metNames(idxDelete) = []; - model.metComps(idxDelete) = []; - model.b(idxDelete) = []; - if isfield(model,'metFormulas') - model.metFormulas(idxDelete) = []; - end - if isfield(model,'unconstrained') - model.unconstrained(idxDelete) = []; - end - if isfield(model,'metMiriams') - model.metMiriams(idxDelete) = []; - end - if isfield(model,'metCharges') - model.metCharges(idxDelete) = []; - end - if isfield(model,'metDeltaG') - model.metDeltaG(idxDelete) = []; - end - if isfield(model,'inchis') - model.inchis(idxDelete) = []; - end - if isfield(model,'metSmiles') - model.metSmiles(idxDelete) = []; - end - if isfield(model,'metFrom') - model.metFrom(idxDelete) = []; - end + repMets = model.mets(repIdx); + model = removeMets(model,idxDelete); + % Indices of the replacement metabolites shift when metabolites with a + % lower index are deleted. + [~, repIdx] = ismember(repMets,model.mets); + repIdx(repIdx==0) = []; end % This could now have created duplicate reactions. Contract model. -model=contractModel(model,[],repIdx); +if isempty(idxDelete) + removedRxns = cell(0,1); + idxDuplRxns = zeros(0,1); +else + originalRxns = model.rxns; + [model, removedRxns] = contractModel(model,[],repIdx); + [~, idxDuplRxns] = ismember(removedRxns,originalRxns); +end end diff --git a/doc/INIT/ftINIT.html b/doc/INIT/ftINIT.html index 6133e894d..8e8a54983 100644 --- a/doc/INIT/ftINIT.html +++ b/doc/INIT/ftINIT.html @@ -270,7 +270,7 @@

SOURCE CODE ^%Get the matching mets 0136 metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); 0137 prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers -0138 (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers +0138 (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' ~= 0); %reversible reactions that are consumers 0139 %convert the production rxns from refModel to minModel 0140 prepData.groupIds; 0141 [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); diff --git a/doc/INIT/mergeLinear.html b/doc/INIT/mergeLinear.html index 16467e80e..5d90456af 100644 --- a/doc/INIT/mergeLinear.html +++ b/doc/INIT/mergeLinear.html @@ -125,157 +125,161 @@

SOURCE CODE ^for i=1:numel(twoNonZero) 0063 involvedRxns=find(reducedModel.S(twoNonZero(i),:)); 0064 %Check that we can have one positive and one negative -0065 pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns)); -0066 neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns)); -0067 -0068 -0069 %Check so that one or both of the reactions haven't been merged -0070 %already -0071 if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1 -0072 %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible, -0073 %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise -0074 if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ... -0075 (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0) -0076 involvedRxns = flip(involvedRxns); -0077 end -0078 -0079 %first make sure the first reaction is producing the metabolite -0080 if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0 -0081 %it is not producing the metabolite - fix that -0082 %first choice: use the second reaction as producer if it is producing -0083 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 -0084 involvedRxns = flip(involvedRxns);%make the second the first -0085 else -0086 %now we know that the second reaction is not producing, so we can safely try to make the first a producer -0087 if reducedModel.rev(involvedRxns(1)) == 1 -0088 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns); -0089 else %ok, finally try to flip the second reaction -0090 if reducedModel.rev(involvedRxns(2)) == 1 -0091 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); -0092 involvedRxns = flip(involvedRxns);%make the second the first -0093 else -0094 error('We should never end up here!'); -0095 end -0096 end -0097 end -0098 end -0099 %Now, make sure the second rxn is a consumer -0100 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 -0101 if reducedModel.rev(involvedRxns(2)) == 1 -0102 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); -0103 else -0104 error('We should never end up here!'); -0105 end -0106 end -0107 -0108 %Calculate how many times the second reaction has to be -0109 %multiplied before being merged with the first -0110 stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2))); +0065 %rev is compared rather than used as-is: the other operand is a +0066 %sparse logical, and MATLAB refuses to combine one with an integer +0067 %array ("Sparse integer array arithmetic operations are not +0068 %supported"). A model built outside RAVEN can carry an integer rev. +0069 pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0); +0070 neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0); +0071 +0072 +0073 %Check so that one or both of the reactions haven't been merged +0074 %already +0075 if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1 +0076 %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible, +0077 %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise +0078 if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ... +0079 (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0) +0080 involvedRxns = flip(involvedRxns); +0081 end +0082 +0083 %first make sure the first reaction is producing the metabolite +0084 if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0 +0085 %it is not producing the metabolite - fix that +0086 %first choice: use the second reaction as producer if it is producing +0087 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 +0088 involvedRxns = flip(involvedRxns);%make the second the first +0089 else +0090 %now we know that the second reaction is not producing, so we can safely try to make the first a producer +0091 if reducedModel.rev(involvedRxns(1)) == 1 +0092 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns); +0093 else %ok, finally try to flip the second reaction +0094 if reducedModel.rev(involvedRxns(2)) == 1 +0095 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); +0096 involvedRxns = flip(involvedRxns);%make the second the first +0097 else +0098 error('We should never end up here!'); +0099 end +0100 end +0101 end +0102 end +0103 %Now, make sure the second rxn is a consumer +0104 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 +0105 if reducedModel.rev(involvedRxns(2)) == 1 +0106 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); +0107 else +0108 error('We should never end up here!'); +0109 end +0110 end 0111 -0112 %Add the second to the first -0113 reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio; -0114 -0115 %Clear the second reaction -0116 reducedModel.S(:,involvedRxns(2))=0; -0117 -0118 %This is to prevent numerical issues. It should be 0 -0119 %already -0120 reducedModel.S(twoNonZero(i),involvedRxns(1))=0; +0112 %Calculate how many times the second reaction has to be +0113 %multiplied before being merged with the first +0114 stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2))); +0115 +0116 %Add the second to the first +0117 reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio; +0118 +0119 %Clear the second reaction +0120 reducedModel.S(:,involvedRxns(2))=0; 0121 -0122 %At this point the second reaction is certain to be deleted -0123 %in a later step and can therefore be ignored -0124 -0125 %Recalculate the bounds for the new reaction. This can be -0126 %problematic since the scale of the bounds may change -0127 %dramatically. Let the most constraining reaction determine -0128 %the new bound -0129 lb1=reducedModel.lb(involvedRxns(1)); -0130 lb2=reducedModel.lb(involvedRxns(2)); -0131 ub1=reducedModel.ub(involvedRxns(1)); -0132 ub2=reducedModel.ub(involvedRxns(2)); -0133 -0134 if lb2~=-inf -0135 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio); -0136 end -0137 if ub2~=inf -0138 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio); -0139 end -0140 -0141 %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev -0142 reducedModel.rev(involvedRxns(1)) = reducedModel.rev(involvedRxns(1))*reducedModel.rev(involvedRxns(2));%this is a way to do an "and" operation with 0 and 1 numbers -0143 -0144 %Then recalculate the objective coefficient. The resulting -0145 %coefficient is the weighted sum of the previous -0146 reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio; +0122 %This is to prevent numerical issues. It should be 0 +0123 %already +0124 reducedModel.S(twoNonZero(i),involvedRxns(1))=0; +0125 +0126 %At this point the second reaction is certain to be deleted +0127 %in a later step and can therefore be ignored +0128 +0129 %Recalculate the bounds for the new reaction. This can be +0130 %problematic since the scale of the bounds may change +0131 %dramatically. Let the most constraining reaction determine +0132 %the new bound +0133 lb1=reducedModel.lb(involvedRxns(1)); +0134 lb2=reducedModel.lb(involvedRxns(2)); +0135 ub1=reducedModel.ub(involvedRxns(1)); +0136 ub2=reducedModel.ub(involvedRxns(2)); +0137 +0138 if lb2~=-inf +0139 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio); +0140 end +0141 if ub2~=inf +0142 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio); +0143 end +0144 +0145 %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev +0146 reducedModel.rev(involvedRxns(1)) = reducedModel.rev(involvedRxns(1))*reducedModel.rev(involvedRxns(2));%this is a way to do an "and" operation with 0 and 1 numbers 0147 -0148 %store which reactions that have been merged -0149 rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1)))); -0150 rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2)))); -0151 grpId = max(groupIds(rxnInd1),groupIds(rxnInd2)); -0152 if grpId == 0 -0153 grpId = nextGroupId; -0154 nextGroupId = nextGroupId + 1; -0155 end -0156 -0157 if groupIds(rxnInd1) ~= grpId -0158 if groupIds(rxnInd1) == 0 -0159 %not merged before, just set the group id -0160 groupIds(rxnInd1) = grpId; -0161 else -0162 %merged before - all rxns with the same group id should be changed -0163 groupIds(groupIds == groupIds(rxnInd1)) = grpId; -0164 end -0165 end -0166 if groupIds(rxnInd2) ~= grpId -0167 if groupIds(rxnInd2) == 0 -0168 %not merged before, just set the group id -0169 groupIds(rxnInd2) = grpId; -0170 else -0171 %merged before - all rxns with the same group id should be changed -0172 groupIds(groupIds == groupIds(rxnInd2)) = grpId; -0173 end -0174 end -0175 -0176 %Iterate again -0177 mergedSome=true; -0178 end -0179 end -0180 -0181 %All possible reactions merged -0182 if mergedSome==false -0183 break; -0184 end -0185 -0186 %Now delete all reactions that involve no metabolites -0187 I=find(sum(reducedModel.S~=0,1)==0); -0188 -0189 %Remove reactions -0190 reducedModel=removeReactions(reducedModel,I); -0191 -0192 %Remove metabolites -0193 notInUse=sum(reducedModel.S~=0,2)==0; -0194 reducedModel=removeMets(reducedModel,notInUse); -0195 end -0196 -0197 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) -0198 model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; -0199 %swap the bounds -0200 ub = model1.ub(rxnInd); -0201 model1.ub(rxnInd) = -model1.lb(rxnInd); -0202 model1.lb(rxnInd) = -ub; -0203 %flip the objective -0204 model1.c(rxnInd) = -model1.c(rxnInd); -0205 -0206 %now take care of the reversedRxns - if this is a group, reverse all of the -0207 %reactions in the group in the reversedRxns index - they will all be reversed at the -0208 %same time since they are the same rxn. -0209 rxnIndices = rxnInd; -0210 if groupIds1(rxnInd) > 0 -0211 rxnIndices = find(groupIds1 == groupIds1(rxnInd)); -0212 end -0213 reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); -0214 end -0215 end +0148 %Then recalculate the objective coefficient. The resulting +0149 %coefficient is the weighted sum of the previous +0150 reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio; +0151 +0152 %store which reactions that have been merged +0153 rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1)))); +0154 rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2)))); +0155 grpId = max(groupIds(rxnInd1),groupIds(rxnInd2)); +0156 if grpId == 0 +0157 grpId = nextGroupId; +0158 nextGroupId = nextGroupId + 1; +0159 end +0160 +0161 if groupIds(rxnInd1) ~= grpId +0162 if groupIds(rxnInd1) == 0 +0163 %not merged before, just set the group id +0164 groupIds(rxnInd1) = grpId; +0165 else +0166 %merged before - all rxns with the same group id should be changed +0167 groupIds(groupIds == groupIds(rxnInd1)) = grpId; +0168 end +0169 end +0170 if groupIds(rxnInd2) ~= grpId +0171 if groupIds(rxnInd2) == 0 +0172 %not merged before, just set the group id +0173 groupIds(rxnInd2) = grpId; +0174 else +0175 %merged before - all rxns with the same group id should be changed +0176 groupIds(groupIds == groupIds(rxnInd2)) = grpId; +0177 end +0178 end +0179 +0180 %Iterate again +0181 mergedSome=true; +0182 end +0183 end +0184 +0185 %All possible reactions merged +0186 if mergedSome==false +0187 break; +0188 end +0189 +0190 %Now delete all reactions that involve no metabolites +0191 I=find(sum(reducedModel.S~=0,1)==0); +0192 +0193 %Remove reactions +0194 reducedModel=removeReactions(reducedModel,I); +0195 +0196 %Remove metabolites +0197 notInUse=sum(reducedModel.S~=0,2)==0; +0198 reducedModel=removeMets(reducedModel,notInUse); +0199 end +0200 +0201 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) +0202 model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; +0203 %swap the bounds +0204 ub = model1.ub(rxnInd); +0205 model1.ub(rxnInd) = -model1.lb(rxnInd); +0206 model1.lb(rxnInd) = -ub; +0207 %flip the objective +0208 model1.c(rxnInd) = -model1.c(rxnInd); +0209 +0210 %now take care of the reversedRxns - if this is a group, reverse all of the +0211 %reactions in the group in the reversedRxns index - they will all be reversed at the +0212 %same time since they are the same rxn. +0213 rxnIndices = rxnInd; +0214 if groupIds1(rxnInd) > 0 +0215 rxnIndices = find(groupIds1 == groupIds1(rxnInd)); +0216 end +0217 reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); +0218 end +0219 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/core/removeMets.html b/doc/core/removeMets.html index 4a3d2d269..2df79cdf2 100644 --- a/doc/core/removeMets.html +++ b/doc/core/removeMets.html @@ -58,7 +58,7 @@

CROSS-REFERENCE INFORMATION ^ <li><a href=convertCharArray convertCharArray
  • getIndexes getIndexes
  • removeReactions removeReactions
    • This function is called by: +
  • mergeCompartments mergeCompartments
  • removeReactions removeReactions
  • replaceMets replaceMets
  • simplifyModel simplifyModel
    • diff --git a/doc/core/replaceMets.html b/doc/core/replaceMets.html index 6f4319513..db6d7f9b3 100644 --- a/doc/core/replaceMets.html +++ b/doc/core/replaceMets.html @@ -48,7 +48,7 @@

    DESCRIPTION ^DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

    This function calls: +
  • contractModel contractModel
  • removeMets removeMets
    • This function is called by: @@ -91,7 +91,7 @@

    SOURCE CODE ^% Output: 0021 % model model structure with selected metabolites replaced 0022 % removedRxns identifiers of duplicate reactions that were removed -0023 % idxDuplRxns index of removedRxns in original model +0023 % idxDuplRxns indices of removedRxns in the input model 0024 % 0025 % Note: This function is useful when the model contains both 'oxygen' and 0026 % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead, @@ -130,97 +130,85 @@

    SOURCE CODE ^'The to-be-replaced metabolite cannot be found in the model.'); 0060 end 0061 -0062 rxnsWithMet = find(model.S(metIdx,:)); -0063 if verbose==true -0064 fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n') -0065 fprintf(strjoin(model.rxns(rxnsWithMet),'\n')) -0066 fprintf('\n') -0067 end -0068 -0069 model.metNames(metIdx) = model.metNames(repIdx(1)); -0070 if isfield(model,'metFormulas') -0071 model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); -0072 end -0073 if isfield(model,'metMiriams') -0074 model.metMiriams(metIdx) = model.metMiriams(repIdx(1)); -0075 end -0076 if isfield(model,'metCharges') -0077 model.metCharges(metIdx) = model.metCharges(repIdx(1)); -0078 end -0079 if isfield(model,'metDeltaG') -0080 model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1)); -0081 end -0082 if isfield(model,'inchis') -0083 model.inchis(metIdx) = model.inchis(repIdx(1)); -0084 end -0085 if isfield(model,'metSmiles') -0086 model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); -0087 end -0088 -0089 idxDelete=[]; -0090 if identifiers -0091 originalStoch = model.S(metIdx,rxnsWithMet); -0092 model.S(repIdx,rxnsWithMet) = originalStoch; -0093 model.S(metIdx,rxnsWithMet) = 0; -0094 idxDelete = metIdx; -0095 else -0096 % Run through replacement metabolites and their compartments. If any of the -0097 % to-be-replaced metabolites is already present (checked by -0098 % metaboliteName[compartment], then the replacement metabolite is kept and -0099 % the to-be-replace metabolite ID deleted. -0100 -0101 % Build list of metaboliteName[compartment] -0102 metCompsN =cellstr(num2str(model.metComps)); -0103 map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps); -0104 metCompsN = map.values(metCompsN); -0105 metCompsN = strcat(lower(model.metNames),'[',metCompsN,']'); -0106 -0107 for i = 1:length(repIdx) -0108 metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN)); -0109 if length(metCompsNidx)>1 -0110 for j = 2:length(metCompsNidx) -0111 model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:); -0112 idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete -0113 end -0114 end -0115 end -0116 end -0117 -0118 if ~isempty(idxDelete) -0119 model.S(idxDelete,:) =[]; -0120 model.mets(idxDelete) = []; -0121 model.metNames(idxDelete) = []; -0122 model.metComps(idxDelete) = []; -0123 model.b(idxDelete) = []; -0124 if isfield(model,'metFormulas') -0125 model.metFormulas(idxDelete) = []; -0126 end -0127 if isfield(model,'unconstrained') -0128 model.unconstrained(idxDelete) = []; -0129 end -0130 if isfield(model,'metMiriams') -0131 model.metMiriams(idxDelete) = []; -0132 end -0133 if isfield(model,'metCharges') -0134 model.metCharges(idxDelete) = []; -0135 end -0136 if isfield(model,'metDeltaG') -0137 model.metDeltaG(idxDelete) = []; -0138 end -0139 if isfield(model,'inchis') -0140 model.inchis(idxDelete) = []; -0141 end -0142 if isfield(model,'metSmiles') -0143 model.metSmiles(idxDelete) = []; -0144 end -0145 if isfield(model,'metFrom') -0146 model.metFrom(idxDelete) = []; -0147 end -0148 end -0149 -0150 % This could now have created duplicate reactions. Contract model. -0151 model=contractModel(model,[],repIdx); -0152 end +0062 [~, rxnsWithMet] = find(model.S(metIdx,:)); +0063 rxnsWithMet = unique(rxnsWithMet); +0064 if verbose==true +0065 fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n') +0066 fprintf('%s\n',model.rxns{rxnsWithMet}) +0067 fprintf('\n') +0068 end +0069 +0070 model.metNames(metIdx) = model.metNames(repIdx(1)); +0071 if isfield(model,'metFormulas') +0072 model.metFormulas(metIdx) = model.metFormulas(repIdx(1)); +0073 end +0074 if isfield(model,'metMiriams') +0075 model.metMiriams(metIdx) = model.metMiriams(repIdx(1)); +0076 end +0077 if isfield(model,'metCharges') +0078 model.metCharges(metIdx) = model.metCharges(repIdx(1)); +0079 end +0080 if isfield(model,'metDeltaG') +0081 model.metDeltaG(metIdx) = model.metDeltaG(repIdx(1)); +0082 end +0083 if isfield(model,'inchis') +0084 model.inchis(metIdx) = model.inchis(repIdx(1)); +0085 end +0086 if isfield(model,'metSmiles') +0087 model.metSmiles(metIdx) = model.metSmiles(repIdx(1)); +0088 end +0089 +0090 idxDelete=[]; +0091 if identifiers +0092 if any(ismember(metIdx,repIdx)) +0093 error('The to-be-replaced metabolite and the replacement metabolite are the same.'); +0094 end +0095 originalStoch = model.S(metIdx,rxnsWithMet); +0096 model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch; +0097 model.S(metIdx,rxnsWithMet) = 0; +0098 idxDelete = metIdx; +0099 else +0100 % Run through replacement metabolites and their compartments. If any of the +0101 % to-be-replaced metabolites is already present (checked by +0102 % metaboliteName[compartment], then the replacement metabolite is kept and +0103 % the to-be-replace metabolite ID deleted. +0104 +0105 % Build list of metaboliteName[compartment] +0106 comps = model.comps(:); +0107 metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']'); +0108 +0109 % Within each metaboliteName[compartment] group, the first replacement +0110 % metabolite is kept and the other members are merged into it. +0111 groups = unique(metCompsN(repIdx),'stable'); +0112 for i = 1:numel(groups) +0113 metCompsNidx = find(strcmp(groups{i}, metCompsN)); +0114 keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1)); +0115 others = metCompsNidx(metCompsNidx ~= keepIdx); +0116 if ~isempty(others) +0117 model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1); +0118 idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete +0119 end +0120 end +0121 end +0122 if ~isempty(idxDelete) +0123 repMets = model.mets(repIdx); +0124 model = removeMets(model,idxDelete); +0125 % Indices of the replacement metabolites shift when metabolites with a +0126 % lower index are deleted. +0127 [~, repIdx] = ismember(repMets,model.mets); +0128 repIdx(repIdx==0) = []; +0129 end +0130 +0131 % This could now have created duplicate reactions. Contract model. +0132 if isempty(idxDelete) +0133 removedRxns = cell(0,1); +0134 idxDuplRxns = zeros(0,1); +0135 else +0136 originalRxns = model.rxns; +0137 [model, removedRxns] = contractModel(model,[],repIdx); +0138 [~, idxDuplRxns] = ismember(removedRxns,originalRxns); +0139 end +0140 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/installation/checkInstallation.html b/doc/installation/checkInstallation.html index 3ee042608..a36271d48 100644 --- a/doc/installation/checkInstallation.html +++ b/doc/installation/checkInstallation.html @@ -381,92 +381,121 @@

    SOURCE CODE ^'hmmerTests.m')"); 0310 res=interpretResults(res); 0311 if res==false -0312 fprintf([' This is essential to run getKEGGModelFromHomology()\n'... -0313 ' when using a FASTA file as input\n']) -0314 end -0315 -0316 if developMode -0317 fprintf('\n=== Development binary executables ===\n'); -0318 fprintf('NOTE: These binaries are only required when using KEGG FTP dump files in getKEGGModelForOrganism\n'); -0319 -0320 fprintf(myStr(' > Checking CD-HIT',40)) -0321 [~,res]=evalc("runtests('cdhitTests.m');"); -0322 interpretResults(res); -0323 -0324 fprintf(myStr(' > Checking MAFFT',40)) -0325 [~,res]=evalc("runtests('mafftTests.m');"); -0326 interpretResults(res); -0327 end -0328 end -0329 -0330 fprintf('\n=== Compatibility ===\n'); -0331 fprintf(myStr(' > Checking function uniqueness',40)) -0332 checkFunctionUniqueness(); -0333 -0334 fprintf('\n*** checkInstallation complete ***\n'); -0335 end -0336 -0337 function res = interpretResults(results) -0338 if results.Failed==0 && results.Incomplete==0 -0339 fprintf('Pass\n'); -0340 res=true; -0341 else -0342 printOrange('Fail\n') -0343 fprintf(' Download/compile the binary and rerun checkInstallation\n'); -0344 res=false; -0345 end -0346 end -0347 -0348 function str = myStr(InputStr,len) -0349 str=InputStr; -0350 lenDiff = len - length(str); -0351 if lenDiff < 0 -0352 warning('String too long'); -0353 else -0354 str = [str blanks(lenDiff)]; -0355 end -0356 end -0357 -0358 function status = makeBinaryExecutable(ravenDir) -0359 % This function is required to run when RAVEN is downloaded as MATLAB -0360 % Add-On, in which case the file permissions are not correctly set -0361 if ispc -0362 status = 0; % No need to run on Windows -0363 return; +0312 fprintf([' HMMER is only required by getKEGGModelForOrganism() when\n'... +0313 ' using a FASTA file as input. If you do not use this function,\n'... +0314 ' this failure can be ignored.\n']) +0315 % On Windows, HMMER runs through Windows Subsystem for Linux +0316 if ispc +0317 [status,~]=system('where wsl'); +0318 if status~=0 +0319 printOrange([' MATLAB cannot find wsl.exe. On Windows, HMMER runs through WSL.\n'... +0320 ' Install it with "wsl --install" in a Command Prompt, then restart\n'... +0321 ' MATLAB (it inherits the system PATH at startup) and rerun checkInstallation\n']) +0322 else +0323 [status,~]=system('wsl -e true'); +0324 if status~=0 +0325 printOrange([' MATLAB finds wsl.exe, but cannot start a Linux distribution\n'... +0326 ' through it. Run "wsl -l -v" in a Command Prompt to check that a\n'... +0327 ' distribution is installed and set up, then restart MATLAB and\n'... +0328 ' rerun checkInstallation\n']) +0329 else +0330 try +0331 [status,~]=system(['wsl "' getWSLpath(fullfile(ravenDir,'software','hmmer','hmmsearch')) '" -h']); +0332 catch +0333 status=1; +0334 end +0335 if status~=0 +0336 printOrange([' WSL works, but the RAVEN HMMER binary cannot be run in it.\n'... +0337 ' Check that the files in software/hmmer are present and that\n'... +0338 ' the drive holding RAVEN is accessible from WSL\n']) +0339 end +0340 end +0341 end +0342 end +0343 end +0344 +0345 if developMode +0346 fprintf('\n=== Development binary executables ===\n'); +0347 fprintf('NOTE: These binaries are only required when using KEGG FTP dump files in getKEGGModelForOrganism\n'); +0348 +0349 fprintf(myStr(' > Checking CD-HIT',40)) +0350 [~,res]=evalc("runtests('cdhitTests.m');"); +0351 interpretResults(res); +0352 +0353 fprintf(myStr(' > Checking MAFFT',40)) +0354 [~,res]=evalc("runtests('mafftTests.m');"); +0355 interpretResults(res); +0356 end +0357 end +0358 +0359 fprintf('\n=== Compatibility ===\n'); +0360 fprintf(myStr(' > Checking function uniqueness',40)) +0361 checkFunctionUniqueness(); +0362 +0363 fprintf('\n*** checkInstallation complete ***\n'); 0364 end -0365 binDir = fullfile(ravenDir,'software'); -0366 -0367 binList = {fullfile(binDir,'blast+','blastp'); fullfile(binDir,'blast+','blastp.mac'); -0368 fullfile(binDir,'blast+','makeblastdb'); fullfile(binDir,'blast+','makeblastdb.mac'); -0369 fullfile(binDir,'cd-hit','cd-hit'); fullfile(binDir,'cd-hit','cd-hit.mac'); -0370 fullfile(binDir,'diamond','diamond'); fullfile(binDir,'diamond','diamond.mac'); -0371 fullfile(binDir,'hmmer','hmmbuild'); fullfile(binDir,'hmmer','hmmbuild.mac'); -0372 fullfile(binDir,'hmmer','hmmsearch'); fullfile(binDir,'hmmer','hmmsearch.mac'); -0373 fullfile(binDir,'GLPKmex','glpkcc.mexa64'); fullfile(binDir,'GLPKmex','glpkcc.mexglx'); fullfile(binDir,'GLPKmex','glpkcc.mexmaci64'); fullfile(binDir,'GLPKmex','glpkcc.mexmaca64'); -0374 fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaca64'); -0375 fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaca64'); -0376 fullfile(binDir,'mafft','mafft-linux64','mafft.bat'); -0377 fullfile(binDir,'mafft','mafft-mac','mafft.bat');}; -0378 -0379 for i=1:numel(binList) -0380 [status,cmdout] = system(['chmod +x "' binList{i} '"']); -0381 if status ~= 0 -0382 warning('Failed to make %s executable: %s ',binList{i},strip(cmdout)) -0383 end +0365 +0366 function res = interpretResults(results) +0367 if results.Failed==0 && results.Incomplete==0 +0368 fprintf('Pass\n'); +0369 res=true; +0370 else +0371 printOrange('Fail\n') +0372 fprintf(' Download/compile the binary and rerun checkInstallation\n'); +0373 res=false; +0374 end +0375 end +0376 +0377 function str = myStr(InputStr,len) +0378 str=InputStr; +0379 lenDiff = len - length(str); +0380 if lenDiff < 0 +0381 warning('String too long'); +0382 else +0383 str = [str blanks(lenDiff)]; 0384 end 0385 end 0386 -0387 function printOrange(stringToPrint) -0388 % printOrange -0389 % Duplicate of RAVEN/core/printOrange is also kept here, as this function -0390 % should be able to run before adding RAVEN to the MATLAB path. -0391 try useDesktop = usejava('desktop'); catch, useDesktop = false; end -0392 if useDesktop -0393 fprintf(['[\b' stringToPrint,']\b']) -0394 else -0395 fprintf(stringToPrint) -0396 end -0397 end +0387 function status = makeBinaryExecutable(ravenDir) +0388 % This function is required to run when RAVEN is downloaded as MATLAB +0389 % Add-On, in which case the file permissions are not correctly set +0390 if ispc +0391 status = 0; % No need to run on Windows +0392 return; +0393 end +0394 binDir = fullfile(ravenDir,'software'); +0395 +0396 binList = {fullfile(binDir,'blast+','blastp'); fullfile(binDir,'blast+','blastp.mac'); +0397 fullfile(binDir,'blast+','makeblastdb'); fullfile(binDir,'blast+','makeblastdb.mac'); +0398 fullfile(binDir,'cd-hit','cd-hit'); fullfile(binDir,'cd-hit','cd-hit.mac'); +0399 fullfile(binDir,'diamond','diamond'); fullfile(binDir,'diamond','diamond.mac'); +0400 fullfile(binDir,'hmmer','hmmbuild'); fullfile(binDir,'hmmer','hmmbuild.mac'); +0401 fullfile(binDir,'hmmer','hmmsearch'); fullfile(binDir,'hmmer','hmmsearch.mac'); +0402 fullfile(binDir,'GLPKmex','glpkcc.mexa64'); fullfile(binDir,'GLPKmex','glpkcc.mexglx'); fullfile(binDir,'GLPKmex','glpkcc.mexmaci64'); fullfile(binDir,'GLPKmex','glpkcc.mexmaca64'); +0403 fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','TranslateSBML_RAVEN.mexmaca64'); +0404 fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexa64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexglx'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaci64'); fullfile(binDir,'libSBML','OutputSBML_RAVEN.mexmaca64'); +0405 fullfile(binDir,'mafft','mafft-linux64','mafft.bat'); +0406 fullfile(binDir,'mafft','mafft-mac','mafft.bat');}; +0407 +0408 for i=1:numel(binList) +0409 [status,cmdout] = system(['chmod +x "' binList{i} '"']); +0410 if status ~= 0 +0411 warning('Failed to make %s executable: %s ',binList{i},strip(cmdout)) +0412 end +0413 end +0414 end +0415 +0416 function printOrange(stringToPrint) +0417 % printOrange +0418 % Duplicate of RAVEN/core/printOrange is also kept here, as this function +0419 % should be able to run before adding RAVEN to the MATLAB path. +0420 try useDesktop = usejava('desktop'); catch, useDesktop = false; end +0421 if useDesktop +0422 fprintf(['[\b' stringToPrint,']\b']) +0423 else +0424 fprintf(stringToPrint) +0425 end +0426 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/importModel.html b/doc/io/importModel.html index dff37de8c..dc124f99e 100644 --- a/doc/io/importModel.html +++ b/doc/io/importModel.html @@ -849,415 +849,426 @@

    SOURCE CODE ^end 0740 end 0741 if isfield(modelSBML,'notes') -0742 startString=strfind(modelSBML.notes,'xhtml">'); -0743 endString=strfind(modelSBML.notes,'</body>'); -0744 if any(startString) && any(endString) -0745 model.annotation.note=modelSBML.notes(startString(1)+7:endString-1); -0746 model.annotation.note=regexprep(model.annotation.note,'<p.*?>|</p.*?>',''); -0747 model.annotation.note=strtrim(model.annotation.note); -0748 if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') -0749 model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel -0750 end -0751 end -0752 end -0753 -0754 if any(~cellfun(@isempty,compartmentOutside)) -0755 model.compOutside=compartmentOutside; -0756 end -0757 -0758 model.rxnNames=reactionNames; -0759 model.metNames=metaboliteNames; -0760 -0761 %Match the compartments for metabolites -0762 [~, J]=ismember(metaboliteCompartments,model.comps); -0763 model.metComps=J; +0742 %Two shapes reach here. RAVEN's own exportModel writes the note as bare +0743 %text inside a <body>; a writer that serialises a notes dictionary emits +0744 %one '<p>key: value</p>' per entry and no <body> at all, carrying the note +0745 %under its own key. Take the labelled entry when there is one, so the +0746 %other entries (version, and whatever else the dictionary held) stay out +0747 %of the note, and fall back to the <body> text otherwise. +0748 note=parseNote(modelSBML.notes,'note'); +0749 if isempty(note) +0750 startString=strfind(modelSBML.notes,'xhtml">'); +0751 endString=strfind(modelSBML.notes,'</body>'); +0752 if any(startString) && any(endString) +0753 note=modelSBML.notes(startString(1)+7:endString-1); +0754 end +0755 end +0756 if ~isempty(note) +0757 model.annotation.note=regexprep(note,'<p.*?>|</p.*?>',''); +0758 model.annotation.note=strtrim(model.annotation.note); +0759 if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') +0760 model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel +0761 end +0762 end +0763 end 0764 -0765 %If any genes have been loaded (only for the new format) -0766 if ~isempty(geneNames) -0767 %In some rare cases geneNames may not necessarily be used in grRules. -0768 %That is true for Yeast 7.6. It's therefore important to change gene -0769 %systematic names to geneIDs in sophisticated way. Gene systematic -0770 %names are not unique, since exactly the same name may be in different -0771 %compartments -0772 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) -0773 geneShortNames=geneNames; -0774 %geneShortNames contain compartments as well, so these are removed -0775 geneShortNames=regexprep(geneShortNames,' \[.+$',''); -0776 %grRules obtained from modifier fields contain geneNames. These are -0777 %changed into geneIDs. grRulesFromModifier is a good way to have -0778 %geneIDs and rxns association when it's important to resolve -0779 %systematic name ambiguities -0780 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0781 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0782 -0783 %Yeast 7.6 contains several metabolites, which were used in gene -0784 %associations. For that reason, the list of species ID is created -0785 %and we then check whether any of them have kegg.genes annotation -0786 %thereby obtaining systematic gene names -0787 geneShortNames=vertcat(geneShortNames,metaboliteNames); -0788 geneIDs=vertcat(geneIDs,metaboliteIDs); -0789 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); -0790 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); -0791 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); -0792 -0793 %Now we retain information for only these entries, which have -0794 %kegg.genes annotation -0795 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); -0796 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); -0797 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); -0798 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); -0799 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); -0800 %Now we reorder geneIDs and geneSystNames by geneSystNames string -0801 %length -0802 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later -0803 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); -0804 geneIDs = geneIDs(Indx); -0805 geneSystNames = geneSystNames(Indx); -0806 for i=1:numel(geneSystNames) -0807 for j=1:numel(grRules) -0808 if strfind(grRules{j},geneSystNames{i}) -0809 if ~isempty(grRules{j}) -0810 if sum(ismember(geneSystNames,geneSystNames{i}))==1 -0811 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); -0812 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 -0813 counter=0; -0814 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); -0815 for k=1:numel(ovrlpIDs) -0816 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) -0817 counter=counter+1; -0818 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); -0819 end -0820 if counter>1 -0821 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; -0822 dispEM(EM); -0823 end -0824 end -0825 end -0826 end -0827 end -0828 end -0829 end -0830 end -0831 model.genes=geneNames; -0832 model.grRules=grRules; -0833 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0834 model.grRules = grRules; -0835 model.rxnGeneMat = rxnGeneMat; -0836 -0837 %Match the compartments for genes -0838 [~, J]=ismember(geneCompartments,model.comps); -0839 model.geneComps=J; -0840 else -0841 if ~all(cellfun(@isempty,grRules)) -0842 %If fbc_geneProduct exists, follow the specified gene order, such -0843 %that matching geneShortNames in function below will work -0844 if isfield(modelSBML,'fbc_geneProduct') -0845 genes={modelSBML.fbc_geneProduct.fbc_id}; -0846 -0847 %Get gene Miriams if they were not retrieved above (this occurs -0848 %when genes are stored as fbc_geneProduct instead of species) -0849 if isempty(geneMiriams) -0850 geneMiriams = cell(numel(genes),1); -0851 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm])) -0852 %If all the SBO terms are identical, don't add them to geneMiriams -0853 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); -0854 end -0855 for i = 1:numel(genes) -0856 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); -0857 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1) -0858 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); -0859 end -0860 end -0861 end -0862 proteins={modelSBML.fbc_geneProduct.fbc_name}; -0863 else -0864 genes=getGenesFromGrRules(grRules); -0865 end -0866 model.genes=genes; -0867 model.grRules=grRules; -0868 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0869 model.grRules = grRules; -0870 model.rxnGeneMat = rxnGeneMat; -0871 end -0872 end -0873 -0874 if all(cellfun(@isempty,geneShortNames)) -0875 if isfield(modelSBML,'fbc_geneProduct') -0876 for i=1:numel(genes) -0877 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) -0878 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; -0879 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) -0880 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; -0881 else -0882 geneShortNames{i,1}=''; -0883 end -0884 end -0885 end -0886 end -0887 -0888 %If any InChIs have been loaded -0889 if any(~cellfun(@isempty,metaboliteInChI)) -0890 model.inchis=metaboliteInChI; -0891 end -0892 -0893 %If any formulas have been loaded -0894 if any(~cellfun(@isempty,metaboliteFormula)) -0895 model.metFormulas=metaboliteFormula; -0896 end -0897 -0898 %If any charges have been loaded -0899 if ~isempty(metaboliteCharges) -0900 model.metCharges=metaboliteCharges; -0901 end -0902 -0903 %If any gene short names have been loaded -0904 if any(~cellfun(@isempty,geneShortNames)) -0905 model.geneShortNames=geneShortNames; -0906 end -0907 -0908 %If any Miriam strings for compartments have been loaded -0909 if any(~cellfun(@isempty,compartmentMiriams)) -0910 model.compMiriams=compartmentMiriams; -0911 end -0912 -0913 %If any Miriam strings for metabolites have been loaded -0914 if any(~cellfun(@isempty,metaboliteMiriams)) -0915 model.metMiriams=metaboliteMiriams; -0916 end -0917 -0918 %If any subsystems have been loaded -0919 if any(~cellfun(@isempty,subsystems)) -0920 model.subSystems=subsystems; -0921 end -0922 if any(rxnComps) -0923 if all(rxnComps) -0924 model.rxnComps=rxnComps; -0925 else -0926 if supressWarnings==false -0927 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; -0928 dispEM(EM,false,model.rxns(rxnComps==0)); -0929 end -0930 end -0931 end -0932 -0933 %If any ec-codes have been loaded -0934 if any(~cellfun(@isempty,eccodes)) -0935 model.eccodes=eccodes; -0936 end -0937 -0938 %If any Miriam strings for reactions have been loaded -0939 if any(~cellfun(@isempty,rxnMiriams)) -0940 model.rxnMiriams=rxnMiriams; -0941 end -0942 -0943 %If any Miriam strings for genes have been loaded -0944 if any(~cellfun(@isempty,geneMiriams)) -0945 model.geneMiriams=geneMiriams; -0946 end -0947 -0948 %If any protein strings have been loaded -0949 if any(~cellfun(@isempty,proteins)) -0950 proteins = reshape(proteins,[],1); -0951 model.proteins=proteins; +0765 if any(~cellfun(@isempty,compartmentOutside)) +0766 model.compOutside=compartmentOutside; +0767 end +0768 +0769 model.rxnNames=reactionNames; +0770 model.metNames=metaboliteNames; +0771 +0772 %Match the compartments for metabolites +0773 [~, J]=ismember(metaboliteCompartments,model.comps); +0774 model.metComps=J; +0775 +0776 %If any genes have been loaded (only for the new format) +0777 if ~isempty(geneNames) +0778 %In some rare cases geneNames may not necessarily be used in grRules. +0779 %That is true for Yeast 7.6. It's therefore important to change gene +0780 %systematic names to geneIDs in sophisticated way. Gene systematic +0781 %names are not unique, since exactly the same name may be in different +0782 %compartments +0783 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) +0784 geneShortNames=geneNames; +0785 %geneShortNames contain compartments as well, so these are removed +0786 geneShortNames=regexprep(geneShortNames,' \[.+$',''); +0787 %grRules obtained from modifier fields contain geneNames. These are +0788 %changed into geneIDs. grRulesFromModifier is a good way to have +0789 %geneIDs and rxns association when it's important to resolve +0790 %systematic name ambiguities +0791 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0792 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0793 +0794 %Yeast 7.6 contains several metabolites, which were used in gene +0795 %associations. For that reason, the list of species ID is created +0796 %and we then check whether any of them have kegg.genes annotation +0797 %thereby obtaining systematic gene names +0798 geneShortNames=vertcat(geneShortNames,metaboliteNames); +0799 geneIDs=vertcat(geneIDs,metaboliteIDs); +0800 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); +0801 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); +0802 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); +0803 +0804 %Now we retain information for only these entries, which have +0805 %kegg.genes annotation +0806 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); +0807 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); +0808 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); +0809 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); +0810 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); +0811 %Now we reorder geneIDs and geneSystNames by geneSystNames string +0812 %length +0813 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later +0814 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); +0815 geneIDs = geneIDs(Indx); +0816 geneSystNames = geneSystNames(Indx); +0817 for i=1:numel(geneSystNames) +0818 for j=1:numel(grRules) +0819 if strfind(grRules{j},geneSystNames{i}) +0820 if ~isempty(grRules{j}) +0821 if sum(ismember(geneSystNames,geneSystNames{i}))==1 +0822 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); +0823 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 +0824 counter=0; +0825 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); +0826 for k=1:numel(ovrlpIDs) +0827 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) +0828 counter=counter+1; +0829 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); +0830 end +0831 if counter>1 +0832 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; +0833 dispEM(EM); +0834 end +0835 end +0836 end +0837 end +0838 end +0839 end +0840 end +0841 end +0842 model.genes=geneNames; +0843 model.grRules=grRules; +0844 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0845 model.grRules = grRules; +0846 model.rxnGeneMat = rxnGeneMat; +0847 +0848 %Match the compartments for genes +0849 [~, J]=ismember(geneCompartments,model.comps); +0850 model.geneComps=J; +0851 else +0852 if ~all(cellfun(@isempty,grRules)) +0853 %If fbc_geneProduct exists, follow the specified gene order, such +0854 %that matching geneShortNames in function below will work +0855 if isfield(modelSBML,'fbc_geneProduct') +0856 genes={modelSBML.fbc_geneProduct.fbc_id}; +0857 +0858 %Get gene Miriams if they were not retrieved above (this occurs +0859 %when genes are stored as fbc_geneProduct instead of species) +0860 if isempty(geneMiriams) +0861 geneMiriams = cell(numel(genes),1); +0862 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && isscalar(unique([modelSBML.fbc_geneProduct.sboTerm])) +0863 %If all the SBO terms are identical, don't add them to geneMiriams +0864 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); +0865 end +0866 for i = 1:numel(genes) +0867 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); +0868 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') && ~(modelSBML.fbc_geneProduct(i).sboTerm==-1) +0869 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); +0870 end +0871 end +0872 end +0873 proteins={modelSBML.fbc_geneProduct.fbc_name}; +0874 else +0875 genes=getGenesFromGrRules(grRules); +0876 end +0877 model.genes=genes; +0878 model.grRules=grRules; +0879 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0880 model.grRules = grRules; +0881 model.rxnGeneMat = rxnGeneMat; +0882 end +0883 end +0884 +0885 if all(cellfun(@isempty,geneShortNames)) +0886 if isfield(modelSBML,'fbc_geneProduct') +0887 for i=1:numel(genes) +0888 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) +0889 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; +0890 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) +0891 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; +0892 else +0893 geneShortNames{i,1}=''; +0894 end +0895 end +0896 end +0897 end +0898 +0899 %If any InChIs have been loaded +0900 if any(~cellfun(@isempty,metaboliteInChI)) +0901 model.inchis=metaboliteInChI; +0902 end +0903 +0904 %If any formulas have been loaded +0905 if any(~cellfun(@isempty,metaboliteFormula)) +0906 model.metFormulas=metaboliteFormula; +0907 end +0908 +0909 %If any charges have been loaded +0910 if ~isempty(metaboliteCharges) +0911 model.metCharges=metaboliteCharges; +0912 end +0913 +0914 %If any gene short names have been loaded +0915 if any(~cellfun(@isempty,geneShortNames)) +0916 model.geneShortNames=geneShortNames; +0917 end +0918 +0919 %If any Miriam strings for compartments have been loaded +0920 if any(~cellfun(@isempty,compartmentMiriams)) +0921 model.compMiriams=compartmentMiriams; +0922 end +0923 +0924 %If any Miriam strings for metabolites have been loaded +0925 if any(~cellfun(@isempty,metaboliteMiriams)) +0926 model.metMiriams=metaboliteMiriams; +0927 end +0928 +0929 %If any subsystems have been loaded +0930 if any(~cellfun(@isempty,subsystems)) +0931 model.subSystems=subsystems; +0932 end +0933 if any(rxnComps) +0934 if all(rxnComps) +0935 model.rxnComps=rxnComps; +0936 else +0937 if supressWarnings==false +0938 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; +0939 dispEM(EM,false,model.rxns(rxnComps==0)); +0940 end +0941 end +0942 end +0943 +0944 %If any ec-codes have been loaded +0945 if any(~cellfun(@isempty,eccodes)) +0946 model.eccodes=eccodes; +0947 end +0948 +0949 %If any Miriam strings for reactions have been loaded +0950 if any(~cellfun(@isempty,rxnMiriams)) +0951 model.rxnMiriams=rxnMiriams; 0952 end 0953 -0954 model.unconstrained=metaboliteUnconstrained; -0955 -0956 %Convert SBML IDs back into their original strings. Here we are using part -0957 %from convertSBMLID, originating from the COBRA Toolbox -0958 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); -0959 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); -0960 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); -0961 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); -0962 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); -0963 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); +0954 %If any Miriam strings for genes have been loaded +0955 if any(~cellfun(@isempty,geneMiriams)) +0956 model.geneMiriams=geneMiriams; +0957 end +0958 +0959 %If any protein strings have been loaded +0960 if any(~cellfun(@isempty,proteins)) +0961 proteins = reshape(proteins,[],1); +0962 model.proteins=proteins; +0963 end 0964 -0965 if removePrefix -0966 [model, hasChanged]=removeIdentifierPrefix(model); -0967 dispEM(['The following fields have prefixes removed from all entries. '... -0968 'If this is undesired, run importModel with removePrefix as false. Example: '... -0969 'importModel(''filename.xml'',[],false);'],false,hasChanged) -0970 end -0971 -0972 %Remove unused fields -0973 if isempty(model.annotation) -0974 model=rmfield(model,'annotation'); -0975 end -0976 if isempty(model.compOutside) -0977 model=rmfield(model,'compOutside'); -0978 end -0979 if isempty(model.compMiriams) -0980 model=rmfield(model,'compMiriams'); +0965 model.unconstrained=metaboliteUnconstrained; +0966 +0967 %Convert SBML IDs back into their original strings. Here we are using part +0968 %from convertSBMLID, originating from the COBRA Toolbox +0969 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); +0970 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); +0971 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); +0972 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); +0973 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); +0974 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); +0975 +0976 if removePrefix +0977 [model, hasChanged]=removeIdentifierPrefix(model); +0978 dispEM(['The following fields have prefixes removed from all entries. '... +0979 'If this is undesired, run importModel with removePrefix as false. Example: '... +0980 'importModel(''filename.xml'',[],false);'],false,hasChanged) 0981 end -0982 if isempty(model.rxnComps) -0983 model=rmfield(model,'rxnComps'); -0984 end -0985 if isempty(model.grRules) -0986 model=rmfield(model,'grRules'); -0987 end -0988 if isempty(model.rxnGeneMat) -0989 model=rmfield(model,'rxnGeneMat'); -0990 end -0991 if isempty(model.subSystems) -0992 model=rmfield(model,'subSystems'); -0993 else -0994 model.subSystems(cellfun(@isempty,subsystems))={{''}}; -0995 % If all subSystems have single entries, then unnest them -0996 if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems)) -0997 model.subSystems = transpose([model.subSystems{:}]); -0998 end -0999 end -1000 if isempty(model.eccodes) -1001 model=rmfield(model,'eccodes'); -1002 end -1003 if isempty(model.rxnMiriams) -1004 model=rmfield(model,'rxnMiriams'); -1005 end -1006 if cellfun(@isempty,model.rxnNotes) -1007 model=rmfield(model,'rxnNotes'); -1008 end -1009 if cellfun(@isempty,model.rxnReferences) -1010 model=rmfield(model,'rxnReferences'); -1011 end -1012 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) -1013 model=rmfield(model,'rxnConfidenceScores'); -1014 end -1015 if isempty(model.genes) -1016 model=rmfield(model,'genes'); -1017 elseif isrow(model.genes) -1018 model.genes=transpose(model.genes); +0982 +0983 %Remove unused fields +0984 if isempty(model.annotation) +0985 model=rmfield(model,'annotation'); +0986 end +0987 if isempty(model.compOutside) +0988 model=rmfield(model,'compOutside'); +0989 end +0990 if isempty(model.compMiriams) +0991 model=rmfield(model,'compMiriams'); +0992 end +0993 if isempty(model.rxnComps) +0994 model=rmfield(model,'rxnComps'); +0995 end +0996 if isempty(model.grRules) +0997 model=rmfield(model,'grRules'); +0998 end +0999 if isempty(model.rxnGeneMat) +1000 model=rmfield(model,'rxnGeneMat'); +1001 end +1002 if isempty(model.subSystems) +1003 model=rmfield(model,'subSystems'); +1004 else +1005 model.subSystems(cellfun(@isempty,subsystems))={{''}}; +1006 % If all subSystems have single entries, then unnest them +1007 if all(cellfun(@(x) iscell(x) && isscalar(x), model.subSystems)) +1008 model.subSystems = transpose([model.subSystems{:}]); +1009 end +1010 end +1011 if isempty(model.eccodes) +1012 model=rmfield(model,'eccodes'); +1013 end +1014 if isempty(model.rxnMiriams) +1015 model=rmfield(model,'rxnMiriams'); +1016 end +1017 if cellfun(@isempty,model.rxnNotes) +1018 model=rmfield(model,'rxnNotes'); 1019 end -1020 if isempty(model.geneComps) -1021 model=rmfield(model,'geneComps'); +1020 if cellfun(@isempty,model.rxnReferences) +1021 model=rmfield(model,'rxnReferences'); 1022 end -1023 if isempty(model.geneMiriams) -1024 model=rmfield(model,'geneMiriams'); +1023 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) +1024 model=rmfield(model,'rxnConfidenceScores'); 1025 end -1026 if isempty(model.geneShortNames) -1027 model=rmfield(model,'geneShortNames'); -1028 end -1029 if isempty(model.proteins) -1030 model=rmfield(model,'proteins'); -1031 end -1032 if isempty(model.inchis) -1033 model=rmfield(model,'inchis'); -1034 end -1035 if isempty(model.metFormulas) -1036 model=rmfield(model,'metFormulas'); -1037 end -1038 if isempty(model.metMiriams) -1039 model=rmfield(model,'metMiriams'); -1040 end -1041 if ~any(model.metCharges) -1042 model=rmfield(model,'metCharges'); -1043 end -1044 -1045 %This just removes the grRules if no genes have been loaded -1046 if ~isfield(model,'genes') && isfield(model,'grRules') -1047 model=rmfield(model,'grRules'); +1026 if isempty(model.genes) +1027 model=rmfield(model,'genes'); +1028 elseif isrow(model.genes) +1029 model.genes=transpose(model.genes); +1030 end +1031 if isempty(model.geneComps) +1032 model=rmfield(model,'geneComps'); +1033 end +1034 if isempty(model.geneMiriams) +1035 model=rmfield(model,'geneMiriams'); +1036 end +1037 if isempty(model.geneShortNames) +1038 model=rmfield(model,'geneShortNames'); +1039 end +1040 if isempty(model.proteins) +1041 model=rmfield(model,'proteins'); +1042 end +1043 if isempty(model.inchis) +1044 model=rmfield(model,'inchis'); +1045 end +1046 if isempty(model.metFormulas) +1047 model=rmfield(model,'metFormulas'); 1048 end -1049 -1050 %Print warnings about bad structure -1051 if supressWarnings==false -1052 checkModelStruct(model,false); -1053 end -1054 -1055 if removeExcMets==true -1056 model=simplifyModel(model); -1057 end -1058 end -1059 -1060 function fieldContent=parseNote(searchString,fieldName) -1061 %The function obtains the particular information from 'notes' field, using -1062 %fieldName as the dummy string -1063 -1064 fieldContent=''; +1049 if isempty(model.metMiriams) +1050 model=rmfield(model,'metMiriams'); +1051 end +1052 if ~any(model.metCharges) +1053 model=rmfield(model,'metCharges'); +1054 end +1055 +1056 %This just removes the grRules if no genes have been loaded +1057 if ~isfield(model,'genes') && isfield(model,'grRules') +1058 model=rmfield(model,'grRules'); +1059 end +1060 +1061 %Print warnings about bad structure +1062 if supressWarnings==false +1063 checkModelStruct(model,false); +1064 end 1065 -1066 if strfind(searchString,fieldName) -1067 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); -1068 targetString=regexprep(targetString,'<p>|</p>',''); -1069 targetString=regexprep(targetString,[fieldName, ':'],''); -1070 for i=1:numel(targetString) -1071 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1072 end -1073 fieldContent=regexprep(fieldContent,'^;|;$',''); -1074 else -1075 fieldContent=''; -1076 end -1077 end -1078 -1079 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) -1080 -1081 fieldContent=''; -1082 -1083 %Removing whitespace characters from the ending strings, which may occur in -1084 %several cases -1085 searchString=regexprep(searchString,'" />','"/>'); -1086 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); -1087 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1088 targetString=regexprep(targetString,startString,''); -1089 targetString=regexprep(targetString,[fieldName midString],''); -1090 -1091 for i=1:numel(targetString) -1092 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1093 end -1094 -1095 fieldContent=regexprep(fieldContent,'^;|;$',''); -1096 end -1097 -1098 function miriamStruct=parseMiriam(searchString) -1099 %Generates miriam structure from annotation field -1100 -1101 %Finding whether miriams are written in the old or the new way -1102 if strfind(searchString,'urn:miriam:') -1103 startString='urn:miriam:'; -1104 elseif strfind(searchString,'http://identifiers.org/') -1105 startString='http://identifiers.org/'; -1106 elseif strfind(searchString,'https://identifiers.org/') -1107 startString='https://identifiers.org/'; -1108 else -1109 miriamStruct=[]; -1110 return; -1111 end -1112 -1113 miriamStruct=[]; -1114 -1115 searchString=regexprep(searchString,'" />','"/>'); -1116 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); -1117 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1118 targetString=regexprep(targetString,startString,''); -1119 -1120 fwdslash = contains(targetString,'/'); -1121 midString = cell(numel(targetString),1); -1122 midString(fwdslash) = {'/'}; -1123 midString(~fwdslash) = {':'}; -1124 -1125 counter=0; -1126 for i=1:numel(targetString) -1127 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) -1128 counter=counter+1; -1129 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); -1130 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); -1131 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); -1132 end -1133 end -1134 end +1066 if removeExcMets==true +1067 model=simplifyModel(model); +1068 end +1069 end +1070 +1071 function fieldContent=parseNote(searchString,fieldName) +1072 %The function obtains the particular information from 'notes' field, using +1073 %fieldName as the dummy string +1074 +1075 fieldContent=''; +1076 +1077 if strfind(searchString,fieldName) +1078 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); +1079 targetString=regexprep(targetString,'<p>|</p>',''); +1080 targetString=regexprep(targetString,[fieldName, ':'],''); +1081 for i=1:numel(targetString) +1082 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1083 end +1084 fieldContent=regexprep(fieldContent,'^;|;$',''); +1085 else +1086 fieldContent=''; +1087 end +1088 end +1089 +1090 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) +1091 +1092 fieldContent=''; +1093 +1094 %Removing whitespace characters from the ending strings, which may occur in +1095 %several cases +1096 searchString=regexprep(searchString,'" />','"/>'); +1097 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); +1098 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1099 targetString=regexprep(targetString,startString,''); +1100 targetString=regexprep(targetString,[fieldName midString],''); +1101 +1102 for i=1:numel(targetString) +1103 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1104 end +1105 +1106 fieldContent=regexprep(fieldContent,'^;|;$',''); +1107 end +1108 +1109 function miriamStruct=parseMiriam(searchString) +1110 %Generates miriam structure from annotation field +1111 +1112 %Finding whether miriams are written in the old or the new way +1113 if strfind(searchString,'urn:miriam:') +1114 startString='urn:miriam:'; +1115 elseif strfind(searchString,'http://identifiers.org/') +1116 startString='http://identifiers.org/'; +1117 elseif strfind(searchString,'https://identifiers.org/') +1118 startString='https://identifiers.org/'; +1119 else +1120 miriamStruct=[]; +1121 return; +1122 end +1123 +1124 miriamStruct=[]; +1125 +1126 searchString=regexprep(searchString,'" />','"/>'); +1127 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); +1128 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1129 targetString=regexprep(targetString,startString,''); +1130 +1131 fwdslash = contains(targetString,'/'); +1132 midString = cell(numel(targetString),1); +1133 midString(fwdslash) = {'/'}; +1134 midString(~fwdslash) = {':'}; 1135 -1136 function miriam = addSBOtoMiriam(miriam,sboTerm) -1137 %Appends SBO term to miriam structure -1138 -1139 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format -1140 if isempty(miriam) -1141 miriam.name = {'sbo'}; -1142 miriam.value = sboTerm; -1143 elseif any(strcmp('sbo',miriam.name)) -1144 currSbo = strcmp('sbo',miriam.name); -1145 miriam.value(currSbo) = sboTerm; -1146 else -1147 miriam.name(end+1) = {'sbo'}; -1148 miriam.value(end+1) = sboTerm; -1149 end -1150 end +1136 counter=0; +1137 for i=1:numel(targetString) +1138 if isempty(regexp(targetString{1,i},'inchi|ec-code', 'once')) +1139 counter=counter+1; +1140 miriamStruct.name{counter,1} = regexprep(targetString{1,i},[midString{i} '.+'],'','once'); +1141 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} midString{i}],'','once'); +1142 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); +1143 end +1144 end +1145 end +1146 +1147 function miriam = addSBOtoMiriam(miriam,sboTerm) +1148 %Appends SBO term to miriam structure +1149 +1150 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format +1151 if isempty(miriam) +1152 miriam.name = {'sbo'}; +1153 miriam.value = sboTerm; +1154 elseif any(strcmp('sbo',miriam.name)) +1155 currSbo = strcmp('sbo',miriam.name); +1156 miriam.value(currSbo) = sboTerm; +1157 else +1158 miriam.name(end+1) = {'sbo'}; +1159 miriam.value(end+1) = sboTerm; +1160 end +1161 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/readYAMLmodel.html b/doc/io/readYAMLmodel.html index d1e03747e..c295b9e19 100644 --- a/doc/io/readYAMLmodel.html +++ b/doc/io/readYAMLmodel.html @@ -205,7 +205,7 @@

    SOURCE CODE ^% import different sections 0150 switch tline_raw -0151 case '- metaData:' +0151 case {'- metaData:', '- metaData: !!omap'} 0152 section = 1; 0153 if verbose 0154 fprintf('\t%d\n', section); @@ -225,7 +225,7 @@

    SOURCE CODE ^end 0169 pos=0; 0170 continue -0171 case '- genes:' +0171 case {'- genes:', '- genes: []'} 0172 section = 4; 0173 if verbose 0174 fprintf('\t%d\n', section); @@ -291,7 +291,7 @@

    SOURCE CODE ^case 'organization' 0236 model.annotation.organization = tline_value; -0237 case 'geckoLight' +0237 case {'geckoLight', 'gecko_light'} 0238 if strcmp(tline_value,'true') 0239 model.ec.geckoLight = true; 0240 end @@ -324,427 +324,466 @@

    SOURCE CODE ^readFieldValue(model, 'inchis', tline_value, pos); 0268 readList=''; miriamKey=''; 0269 case 'smiles' -0270 model = readFieldValue(model, 'metSmiles', tline_value, pos); -0271 readList=''; miriamKey=''; -0272 case 'deltaG' -0273 model = readFieldValue(model, 'metDeltaG', tline_value, pos); -0274 readList=''; miriamKey=''; -0275 case 'metFrom' -0276 model = readFieldValue(model, 'metFrom', tline_value, pos); -0277 readList=''; miriamKey=''; -0278 case 'annotation' -0279 readList = 'annotation'; -0280 otherwise -0281 switch readList -0282 case 'annotation' -0283 [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); -0284 otherwise -0285 error(['Unknown entry in yaml file: ' tline_raw]) -0286 end -0287 end; continue -0288 end -0289 -0290 % import reactions: -0291 if section == 3 -0292 switch tline_key -0293 case 'id' -0294 pos = pos + 1; -0295 model = readFieldValue(model, 'rxns', tline_value,pos); -0296 readList=''; miriamKey=''; -0297 case 'name' -0298 model = readFieldValue(model, 'rxnNames', tline_value, pos); -0299 readList=''; miriamKey=''; -0300 case 'lower_bound' -0301 model.lb(pos,1) = {tline_value}; -0302 readList=''; miriamKey=''; -0303 case 'upper_bound' -0304 model.ub(pos,1) = {tline_value}; -0305 readList=''; miriamKey=''; -0306 case 'rev' -0307 model.rev(pos,1) = {tline_value}; -0308 readList=''; miriamKey=''; -0309 case 'gene_reaction_rule' -0310 model = readFieldValue(model, 'grRules', tline_value, pos); -0311 readList=''; miriamKey=''; -0312 case 'rxnNotes' -0313 model = readFieldValue(model, 'rxnNotes', tline_value, pos); -0314 readList=''; miriamKey=''; -0315 case 'rxnFrom' -0316 model = readFieldValue(model, 'rxnFrom', tline_value, pos); -0317 readList=''; miriamKey=''; -0318 case 'deltaG' -0319 model = readFieldValue(model, 'rxnDeltaG', tline_value, pos); -0320 readList=''; miriamKey=''; -0321 case 'objective_coefficient' -0322 model.c(pos,1) = 1; -0323 readList=''; miriamKey=''; -0324 case 'references' -0325 model = readFieldValue(model, 'rxnReferences', tline_value, pos); -0326 readList=''; miriamKey=''; -0327 case 'confidence_score' -0328 model = readFieldValue(model, 'rxnConfidenceScores', tline_value, pos); -0329 readList=''; miriamKey=''; -0330 case 'eccodes' -0331 if isempty(tline_value) -0332 readList = 'eccodes'; -0333 else -0334 eccodes(ecCodeNo,1:2)={pos,tline_value}; -0335 ecCodeNo=ecCodeNo+1; -0336 end -0337 case 'subsystem' -0338 if isempty(tline_value) -0339 readList = 'subsystem'; -0340 else -0341 subSystems(subSysNo,1:2)={pos,tline_value}; -0342 subSysNo=subSysNo+1; -0343 end -0344 case 'metabolites' -0345 readList = 'equation'; -0346 case 'annotation' -0347 readList = 'annotation'; -0348 -0349 otherwise -0350 switch readList -0351 case 'eccodes' -0352 eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; -0353 ecCodeNo=ecCodeNo+1; -0354 case 'subsystem' -0355 subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; -0356 subSysNo=subSysNo+1; -0357 case 'annotation' -0358 [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); -0359 case 'equation' -0360 coeff = sscanf(tline_value,'%f'); -0361 equations(equatiNo,1:3)={pos,tline_key,coeff}; -0362 equatiNo=equatiNo+1; -0363 otherwise -0364 error(['Unknown entry in yaml file: ' tline_raw]) -0365 end -0366 end; continue -0367 end -0368 -0369 % import genes: -0370 if section == 4 -0371 switch tline_key -0372 case 'id' -0373 pos = pos + 1; -0374 model = readFieldValue(model, 'genes', tline_value, pos); -0375 readList = ''; -0376 miriamKey = ''; -0377 case 'name' -0378 model = readFieldValue(model, 'geneShortNames', tline_value, pos); -0379 case 'protein' -0380 model = readFieldValue(model, 'proteins', tline_value, pos); -0381 case 'annotation' -0382 readList = 'annotation'; -0383 otherwise -0384 switch readList -0385 case 'annotation' -0386 [geneMiriams, miriamKey,genMirNo] = gatherAnnotation(pos,geneMiriams,tline_key,tline_value,miriamKey,genMirNo); -0387 otherwise -0388 error(['Unknown entry in yaml file: ' tline_raw]) -0389 end -0390 end; continue -0391 end -0392 -0393 % import compartments: -0394 if section == 5 -0395 model.comps(end+1,1) = {tline_key}; -0396 model.compNames(end+1,1) = {tline_value}; -0397 end -0398 -0399 % import ec reaction info -0400 if section == 6 -0401 switch tline_key -0402 case 'id' -0403 pos = pos + 1; -0404 model.ec = readFieldValue(model.ec, 'rxns', tline_value, pos); -0405 readList=''; -0406 case 'kcat' -0407 model.ec = readFieldValue(model.ec, 'kcat', tline_value, pos); -0408 readList=''; -0409 case 'source' -0410 model.ec = readFieldValue(model.ec, 'source', tline_value, pos); -0411 readList=''; -0412 case 'notes' -0413 model.ec = readFieldValue(model.ec, 'notes', tline_value, pos); -0414 readList=''; -0415 case 'eccodes' -0416 if isempty(tline_value) -0417 readList = 'eccodes'; -0418 else -0419 ecGecko(ecGeckoNo,1:2)={pos,tline_value}; -0420 ecGeckoNo=ecGeckoNo+1; -0421 end -0422 case 'enzymes' -0423 readList = 'enzStoich'; -0424 otherwise -0425 switch readList -0426 case 'eccodes' -0427 ecGecko(ecGeckoNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; -0428 ecGeckoNo=ecGeckoNo+1; -0429 case 'enzStoich' -0430 coeff = sscanf(tline_value,'%f'); -0431 enzStoich(enzStoichNo,1:3)={pos,tline_key,coeff}; -0432 enzStoichNo=enzStoichNo+1; -0433 otherwise -0434 error(['Unknown entry in yaml file: ' tline_raw]) -0435 end -0436 end; continue -0437 end -0438 -0439 % import ec enzyme info -0440 if section == 7 -0441 switch tline_key -0442 case 'genes' -0443 pos = pos + 1; -0444 model.ec = readFieldValue(model.ec, 'genes', tline_value, pos); -0445 case 'enzymes' -0446 model.ec = readFieldValue(model.ec, 'enzymes', tline_value, pos); -0447 case 'mw' -0448 model.ec = readFieldValue(model.ec, 'mw', tline_value, pos); -0449 case 'sequence' -0450 model.ec = readFieldValue(model.ec, 'sequence', tline_value, pos); -0451 case 'concs' -0452 model.ec = readFieldValue(model.ec, 'concs', tline_value, pos); -0453 otherwise -0454 error(['Unknown entry in yaml file: ' tline_raw]) -0455 end; continue -0456 end -0457 end -0458 -0459 %Parse annotations -0460 if ~isempty(metMiriams) -0461 locs = cell2mat(metMiriams(:,1)); -0462 for i=unique(locs)' -0463 model.metMiriams{i,1}.name=metMiriams(locs==i,2); -0464 model.metMiriams{i,1}.value=metMiriams(locs==i,3); -0465 end -0466 end -0467 if ~isempty(rxnMiriams) -0468 locs = cell2mat(rxnMiriams(:,1)); -0469 for i=unique(locs)' -0470 model.rxnMiriams{i,1}.name=rxnMiriams(locs==i,2); -0471 model.rxnMiriams{i,1}.value=rxnMiriams(locs==i,3); -0472 end -0473 end -0474 if ~isempty(geneMiriams) -0475 locs = cell2mat(geneMiriams(:,1)); -0476 for i=unique(locs)' -0477 model.geneMiriams{i,1}.name=geneMiriams(locs==i,2); -0478 model.geneMiriams{i,1}.value=geneMiriams(locs==i,3); -0479 end -0480 end -0481 -0482 %Parse subSystems -0483 if ~isempty(subSystems) -0484 locs = cell2mat(subSystems(:,1)); -0485 for i=unique(locs)' -0486 model.subSystems{i,1}=subSystems(locs==i,2); -0487 end -0488 end -0489 -0490 %Parse ec-codes -0491 if ~isempty(eccodes) -0492 locs = cell2mat(eccodes(:,1)); -0493 for i=unique(locs)' -0494 eccodesCat=strjoin(eccodes(locs==i,2),';'); -0495 model.eccodes{i,1}=eccodesCat; -0496 end -0497 emptyEc=cellfun('isempty',model.eccodes); -0498 model.eccodes(emptyEc)={''}; -0499 end -0500 -0501 % follow-up data processing -0502 if verbose -0503 fprintf('\nimporting completed\nfollow-up processing...'); -0504 end -0505 [~, model.metComps] = ismember(model.metComps, model.comps); -0506 [~, model.geneComps] = ismember(model.geneComps, model.comps); -0507 [~, model.rxnComps] = ismember(model.rxnComps, model.comps); -0508 -0509 % Fill S-matrix -0510 rxnIdx = cellfun('isempty', equations(:,1)); -0511 equations(rxnIdx,:) = ''; -0512 rxnIdx = cell2mat(equations(:,1)); -0513 [~,metIdx] = ismember(equations(:,2),model.mets); -0514 coeffs = cell2mat(equations(:,3)); -0515 model.S=sparse(max(metIdx),max(rxnIdx)); -0516 linearIndices = sub2ind([max(metIdx), max(rxnIdx)],metIdx,rxnIdx); -0517 model.S(linearIndices) = coeffs; -0518 -0519 % Convert strings to numeric -0520 model.metCharges = str2double(model.metCharges); -0521 model.lb = str2double(model.lb); -0522 model.ub = str2double(model.ub); -0523 model.rxnConfidenceScores = str2double(model.rxnConfidenceScores); -0524 model.b = zeros(length(model.mets),1); -0525 model.metDeltaG = str2double(model.metDeltaG); -0526 model.rxnDeltaG = str2double(model.rxnDeltaG); -0527 -0528 % Fill some other fields -0529 model.annotation.defaultLB = min(model.lb); -0530 model.annotation.defaultUB = max(model.ub); -0531 if numel(model.lb)<numel(model.rxns) %No LB reported = min -0532 model.lb(end+1:numel(model.rxns)-numel(model.lb),1) = double(model.annotation.defaultLB); -0533 end -0534 if numel(model.ub)<numel(model.rxns) %No UB reported = max -0535 model.ub(end+1:numel(model.rxns)-numel(model.ub),1) = double(model.annotation.defaultUB); -0536 end -0537 if ~all(cellfun('isempty',model.rev)) -0538 model.rev = str2double(model.rev); -0539 else -0540 model.rev = []; -0541 end -0542 if numel(model.rev)<numel(model.rxns) %No rev reported, assume from LB and UB -0543 model.rev(end+1:numel(model.rxns)-numel(model.rev),1) = double(model.lb<0 & model.ub>0); -0544 end -0545 -0546 % Remove empty fields, otherwise fill to correct length -0547 % Reactions -0548 for i={'rxnNames','grRules','eccodes','rxnNotes','rxnReferences',... -0549 'rxnFrom','subSystems','rxnMiriams'} % Empty strings -0550 model = emptyOrFill(model,i{1},{''},'rxns'); -0551 end -0552 for i={'c'} % Zeros -0553 model = emptyOrFill(model,i{1},0,'rxns',true); -0554 end -0555 for i={'rxnConfidenceScores','rxnDeltaG'} % NaNs -0556 model = emptyOrFill(model,i{1},NaN,'rxns'); -0557 end -0558 for i={'rxnComps'} % Ones, assume first compartment -0559 model = emptyOrFill(model,i{1},1,'rxns'); -0560 end -0561 % Metabolites -0562 for i={'metNames','inchis','metFormulas','metMiriams','metFrom','metSmiles','metNotes'} % Empty strings -0563 model = emptyOrFill(model,i{1},{''},'mets'); -0564 end -0565 for i={'metCharges','unconstrained'} % Zeros -0566 model = emptyOrFill(model,i{1},0,'mets'); -0567 end -0568 for i={'metDeltaG'} % % NaNs -0569 model = emptyOrFill(model,i{1},NaN,'mets'); -0570 end -0571 for i={'metComps'} % Ones, assume first compartment -0572 model = emptyOrFill(model,i{1},1,'mets'); -0573 end -0574 % Genes -0575 for i={'geneMiriams','geneShortNames','proteins'} % Empty strings -0576 model = emptyOrFill(model,i{1},{''},'genes'); -0577 end -0578 for i={'geneComps'} % Ones, assume first compartment -0579 model = emptyOrFill(model,i{1},1,'genes'); +0270 if isempty(tline_value) +0271 readList = 'smilesInAnnotation'; +0272 else +0273 model = readFieldValue(model, 'metSmiles', tline_value, pos); +0274 readList=''; miriamKey=''; +0275 end +0276 case 'deltaG' +0277 model = readFieldValue(model, 'metDeltaG', tline_value, pos); +0278 readList=''; miriamKey=''; +0279 case 'metFrom' +0280 model = readFieldValue(model, 'metFrom', tline_value, pos); +0281 readList=''; miriamKey=''; +0282 case 'annotation' +0283 readList = 'annotation'; +0284 otherwise +0285 switch readList +0286 case 'smilesInAnnotation' +0287 %As for ec-code on a reaction: a bare list item +0288 %continues the list, anything else ends it. +0289 if isempty(tline_key) +0290 model = readFieldValue(model, 'metSmiles', ... +0291 regexprep(tline_value,'^ +- "?(.*)"?$','$1'), pos); +0292 else +0293 readList = 'annotation'; +0294 [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); +0295 end +0296 case 'annotation' +0297 [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); +0298 otherwise +0299 error(['Unknown entry in yaml file: ' tline_raw]) +0300 end +0301 end; continue +0302 end +0303 +0304 % import reactions: +0305 if section == 3 +0306 switch tline_key +0307 case 'id' +0308 pos = pos + 1; +0309 model = readFieldValue(model, 'rxns', tline_value,pos); +0310 readList=''; miriamKey=''; +0311 case 'name' +0312 model = readFieldValue(model, 'rxnNames', tline_value, pos); +0313 readList=''; miriamKey=''; +0314 case 'lower_bound' +0315 model.lb(pos,1) = {tline_value}; +0316 readList=''; miriamKey=''; +0317 case 'upper_bound' +0318 model.ub(pos,1) = {tline_value}; +0319 readList=''; miriamKey=''; +0320 case 'rev' +0321 model.rev(pos,1) = {tline_value}; +0322 readList=''; miriamKey=''; +0323 case 'gene_reaction_rule' +0324 model = readFieldValue(model, 'grRules', tline_value, pos); +0325 readList=''; miriamKey=''; +0326 case {'rxnNotes', 'notes'} +0327 model = readFieldValue(model, 'rxnNotes', tline_value, pos); +0328 readList=''; miriamKey=''; +0329 case 'rxnFrom' +0330 model = readFieldValue(model, 'rxnFrom', tline_value, pos); +0331 readList=''; miriamKey=''; +0332 case 'deltaG' +0333 model = readFieldValue(model, 'rxnDeltaG', tline_value, pos); +0334 readList=''; miriamKey=''; +0335 case 'objective_coefficient' +0336 model.c(pos,1) = 1; +0337 readList=''; miriamKey=''; +0338 case 'references' +0339 model = readFieldValue(model, 'rxnReferences', tline_value, pos); +0340 readList=''; miriamKey=''; +0341 case 'confidence_score' +0342 model = readFieldValue(model, 'rxnConfidenceScores', tline_value, pos); +0343 readList=''; miriamKey=''; +0344 case 'eccodes' +0345 if isempty(tline_value) +0346 readList = 'eccodes'; +0347 else +0348 eccodes(ecCodeNo,1:2)={pos,tline_value}; +0349 ecCodeNo=ecCodeNo+1; +0350 end +0351 case 'subsystem' +0352 if isempty(tline_value) +0353 readList = 'subsystem'; +0354 else +0355 subSystems(subSysNo,1:2)={pos,tline_value}; +0356 subSysNo=subSysNo+1; +0357 end +0358 case 'metabolites' +0359 readList = 'equation'; +0360 case 'ec-code' +0361 %EC numbers emitted inside the annotation block rather than +0362 %under the legacy top-level 'eccodes' key, which is still +0363 %accepted above. A block list follows when the value is +0364 %empty; it is read under its own state so that gathering can +0365 %resume for the rest of the annotation block once the list +0366 %ends, instead of swallowing every later entry as an EC code. +0367 if isempty(tline_value) +0368 readList = 'eccodesInAnnotation'; +0369 else +0370 eccodes(ecCodeNo,1:2)={pos,tline_value}; +0371 ecCodeNo=ecCodeNo+1; +0372 end +0373 case 'annotation' +0374 readList = 'annotation'; +0375 +0376 otherwise +0377 switch readList +0378 case 'eccodes' +0379 eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0380 ecCodeNo=ecCodeNo+1; +0381 case 'eccodesInAnnotation' +0382 %A bare list item ("- 1.1.2.4") has no colon, so +0383 %tline_key is empty and the list continues. Anything +0384 %else is the annotation block's next entry, so the +0385 %list has ended. +0386 if isempty(tline_key) +0387 eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0388 ecCodeNo=ecCodeNo+1; +0389 else +0390 readList = 'annotation'; +0391 [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); +0392 end +0393 case 'subsystem' +0394 subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0395 subSysNo=subSysNo+1; +0396 case 'annotation' +0397 [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); +0398 case 'equation' +0399 coeff = sscanf(tline_value,'%f'); +0400 equations(equatiNo,1:3)={pos,tline_key,coeff}; +0401 equatiNo=equatiNo+1; +0402 otherwise +0403 error(['Unknown entry in yaml file: ' tline_raw]) +0404 end +0405 end; continue +0406 end +0407 +0408 % import genes: +0409 if section == 4 +0410 switch tline_key +0411 case 'id' +0412 pos = pos + 1; +0413 model = readFieldValue(model, 'genes', tline_value, pos); +0414 readList = ''; +0415 miriamKey = ''; +0416 case 'name' +0417 model = readFieldValue(model, 'geneShortNames', tline_value, pos); +0418 case 'protein' +0419 model = readFieldValue(model, 'proteins', tline_value, pos); +0420 case 'annotation' +0421 readList = 'annotation'; +0422 otherwise +0423 switch readList +0424 case 'annotation' +0425 [geneMiriams, miriamKey,genMirNo] = gatherAnnotation(pos,geneMiriams,tline_key,tline_value,miriamKey,genMirNo); +0426 otherwise +0427 error(['Unknown entry in yaml file: ' tline_raw]) +0428 end +0429 end; continue +0430 end +0431 +0432 % import compartments: +0433 if section == 5 +0434 model.comps(end+1,1) = {tline_key}; +0435 model.compNames(end+1,1) = {tline_value}; +0436 end +0437 +0438 % import ec reaction info +0439 if section == 6 +0440 switch tline_key +0441 case 'id' +0442 pos = pos + 1; +0443 model.ec = readFieldValue(model.ec, 'rxns', tline_value, pos); +0444 readList=''; +0445 case 'kcat' +0446 model.ec = readFieldValue(model.ec, 'kcat', tline_value, pos); +0447 readList=''; +0448 case 'source' +0449 model.ec = readFieldValue(model.ec, 'source', tline_value, pos); +0450 readList=''; +0451 case 'notes' +0452 model.ec = readFieldValue(model.ec, 'notes', tline_value, pos); +0453 readList=''; +0454 case 'eccodes' +0455 if isempty(tline_value) +0456 readList = 'eccodes'; +0457 else +0458 ecGecko(ecGeckoNo,1:2)={pos,tline_value}; +0459 ecGeckoNo=ecGeckoNo+1; +0460 end +0461 case 'enzymes' +0462 readList = 'enzStoich'; +0463 otherwise +0464 switch readList +0465 case 'eccodes' +0466 ecGecko(ecGeckoNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; +0467 ecGeckoNo=ecGeckoNo+1; +0468 case 'enzStoich' +0469 coeff = sscanf(tline_value,'%f'); +0470 enzStoich(enzStoichNo,1:3)={pos,tline_key,coeff}; +0471 enzStoichNo=enzStoichNo+1; +0472 otherwise +0473 error(['Unknown entry in yaml file: ' tline_raw]) +0474 end +0475 end; continue +0476 end +0477 +0478 % import ec enzyme info +0479 if section == 7 +0480 switch tline_key +0481 case 'genes' +0482 pos = pos + 1; +0483 model.ec = readFieldValue(model.ec, 'genes', tline_value, pos); +0484 case 'enzymes' +0485 model.ec = readFieldValue(model.ec, 'enzymes', tline_value, pos); +0486 case 'mw' +0487 model.ec = readFieldValue(model.ec, 'mw', tline_value, pos); +0488 case 'sequence' +0489 model.ec = readFieldValue(model.ec, 'sequence', tline_value, pos); +0490 case 'concs' +0491 model.ec = readFieldValue(model.ec, 'concs', tline_value, pos); +0492 otherwise +0493 error(['Unknown entry in yaml file: ' tline_raw]) +0494 end; continue +0495 end +0496 end +0497 +0498 %Parse annotations +0499 if ~isempty(metMiriams) +0500 locs = cell2mat(metMiriams(:,1)); +0501 for i=unique(locs)' +0502 model.metMiriams{i,1}.name=metMiriams(locs==i,2); +0503 model.metMiriams{i,1}.value=metMiriams(locs==i,3); +0504 end +0505 end +0506 if ~isempty(rxnMiriams) +0507 locs = cell2mat(rxnMiriams(:,1)); +0508 for i=unique(locs)' +0509 model.rxnMiriams{i,1}.name=rxnMiriams(locs==i,2); +0510 model.rxnMiriams{i,1}.value=rxnMiriams(locs==i,3); +0511 end +0512 end +0513 if ~isempty(geneMiriams) +0514 locs = cell2mat(geneMiriams(:,1)); +0515 for i=unique(locs)' +0516 model.geneMiriams{i,1}.name=geneMiriams(locs==i,2); +0517 model.geneMiriams{i,1}.value=geneMiriams(locs==i,3); +0518 end +0519 end +0520 +0521 %Parse subSystems +0522 if ~isempty(subSystems) +0523 locs = cell2mat(subSystems(:,1)); +0524 for i=unique(locs)' +0525 model.subSystems{i,1}=subSystems(locs==i,2); +0526 end +0527 end +0528 +0529 %Parse ec-codes +0530 if ~isempty(eccodes) +0531 locs = cell2mat(eccodes(:,1)); +0532 for i=unique(locs)' +0533 eccodesCat=strjoin(eccodes(locs==i,2),';'); +0534 model.eccodes{i,1}=eccodesCat; +0535 end +0536 emptyEc=cellfun('isempty',model.eccodes); +0537 model.eccodes(emptyEc)={''}; +0538 end +0539 +0540 % follow-up data processing +0541 if verbose +0542 fprintf('\nimporting completed\nfollow-up processing...'); +0543 end +0544 [~, model.metComps] = ismember(model.metComps, model.comps); +0545 [~, model.geneComps] = ismember(model.geneComps, model.comps); +0546 [~, model.rxnComps] = ismember(model.rxnComps, model.comps); +0547 +0548 % Fill S-matrix +0549 rxnIdx = cellfun('isempty', equations(:,1)); +0550 equations(rxnIdx,:) = ''; +0551 rxnIdx = cell2mat(equations(:,1)); +0552 [~,metIdx] = ismember(equations(:,2),model.mets); +0553 coeffs = cell2mat(equations(:,3)); +0554 model.S=sparse(max(metIdx),max(rxnIdx)); +0555 linearIndices = sub2ind([max(metIdx), max(rxnIdx)],metIdx,rxnIdx); +0556 model.S(linearIndices) = coeffs; +0557 +0558 % Convert strings to numeric +0559 model.metCharges = str2double(model.metCharges); +0560 model.lb = str2double(model.lb); +0561 model.ub = str2double(model.ub); +0562 model.rxnConfidenceScores = str2double(model.rxnConfidenceScores); +0563 model.b = zeros(length(model.mets),1); +0564 model.metDeltaG = str2double(model.metDeltaG); +0565 model.rxnDeltaG = str2double(model.rxnDeltaG); +0566 +0567 % Fill some other fields +0568 model.annotation.defaultLB = min(model.lb); +0569 model.annotation.defaultUB = max(model.ub); +0570 if numel(model.lb)<numel(model.rxns) %No LB reported = min +0571 model.lb(end+1:numel(model.rxns)-numel(model.lb),1) = double(model.annotation.defaultLB); +0572 end +0573 if numel(model.ub)<numel(model.rxns) %No UB reported = max +0574 model.ub(end+1:numel(model.rxns)-numel(model.ub),1) = double(model.annotation.defaultUB); +0575 end +0576 if ~all(cellfun('isempty',model.rev)) +0577 model.rev = str2double(model.rev); +0578 else +0579 model.rev = []; 0580 end -0581 % Comps -0582 for i={'compNames'} % Empty strings -0583 model = emptyOrFill(model,i{1},{''},'comps'); -0584 end -0585 for i={'compOutside'} % First comp -0586 model = emptyOrFill(model,i{1},model.comps{1},'comps'); -0587 end -0588 % Single fields are kept, even if empty -0589 % for i={'description','name','version','date','annotation'} -0590 % if isempty(model.(i{1})) -0591 % model = rmfield(model,i{1}); -0592 % end -0593 % end -0594 -0595 if isfield(model,'subSystems') -0596 % If all entries are 1x1, then flatten -0597 if all(cellfun(@(x) numel(x) <= 1, model.subSystems)) -0598 model.subSystems = transpose([model.subSystems{:}]); -0599 end -0600 end -0601 -0602 % Make rxnGeneMat fields and map to the existing model.genes field -0603 if isfield(model,'grRules') -0604 [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules); -0605 model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes)); -0606 [~,geneOrder] = ismember(genes,model.genes); -0607 if any(geneOrder == 0) -0608 error(['The grRules includes the following gene(s), that are not in '... -0609 'the list of model genes: ', genes{~geneOrder}]) -0610 end -0611 model.rxnGeneMat(:,geneOrder) = rxnGeneMat; +0581 if numel(model.rev)<numel(model.rxns) %No rev reported, assume from LB and UB +0582 model.rev(end+1:numel(model.rxns)-numel(model.rev),1) = double(model.lb<0 & model.ub>0); +0583 end +0584 +0585 % Remove empty fields, otherwise fill to correct length +0586 % Reactions +0587 for i={'rxnNames','grRules','eccodes','rxnNotes','rxnReferences',... +0588 'rxnFrom','subSystems','rxnMiriams'} % Empty strings +0589 model = emptyOrFill(model,i{1},{''},'rxns'); +0590 end +0591 for i={'c'} % Zeros +0592 model = emptyOrFill(model,i{1},0,'rxns',true); +0593 end +0594 for i={'rxnConfidenceScores','rxnDeltaG'} % NaNs +0595 model = emptyOrFill(model,i{1},NaN,'rxns'); +0596 end +0597 for i={'rxnComps'} % Ones, assume first compartment +0598 model = emptyOrFill(model,i{1},1,'rxns'); +0599 end +0600 % Metabolites +0601 for i={'metNames','inchis','metFormulas','metMiriams','metFrom','metSmiles','metNotes'} % Empty strings +0602 model = emptyOrFill(model,i{1},{''},'mets'); +0603 end +0604 for i={'metCharges','unconstrained'} % Zeros +0605 model = emptyOrFill(model,i{1},0,'mets'); +0606 end +0607 for i={'metDeltaG'} % % NaNs +0608 model = emptyOrFill(model,i{1},NaN,'mets'); +0609 end +0610 for i={'metComps'} % Ones, assume first compartment +0611 model = emptyOrFill(model,i{1},1,'mets'); 0612 end -0613 -0614 % Finalize GECKO model -0615 if isGECKO -0616 % Fill in empty fields and empty entries -0617 for i={'kcat','source','notes','eccodes'} % Even keep empty -0618 model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true); -0619 end -0620 for i={'enzymes','mw','sequence'} -0621 model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true); -0622 end -0623 model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true); -0624 model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true); -0625 % Change string to double -0626 for i={'kcat','mw','concs'} -0627 if isfield(model.ec,i{1}) -0628 model.ec.(i{1}) = str2double(model.ec.(i{1})); -0629 end -0630 end -0631 % Fill rxnEnzMat -0632 rxnIdx = cellfun('isempty', enzStoich(:,1)); -0633 enzStoich(rxnIdx,:) = ''; -0634 rxnIdx = cell2mat(enzStoich(:,1)); -0635 [~,enzIdx] = ismember(enzStoich(:,2),model.ec.enzymes); -0636 coeffs = cell2mat(enzStoich(:,3)); -0637 model.ec.rxnEnzMat = zeros(numel(model.ec.rxns), numel(model.ec.genes)); -0638 linearIndices = sub2ind([numel(model.ec.rxns), numel(model.ec.genes)], rxnIdx, enzIdx); -0639 model.ec.rxnEnzMat(linearIndices) = coeffs; -0640 %Parse ec-codes -0641 if ~isempty(ecGecko) -0642 locs = cell2mat(ecGecko(:,1)); -0643 for i=unique(locs)' -0644 ecGeckoCat=strjoin(ecGecko(locs==i,2),';'); -0645 model.ec.eccodes{i,1}=ecGeckoCat; -0646 end -0647 emptyEc=cellfun('isempty',model.ec.eccodes); -0648 model.ec.eccodes(emptyEc)={''}; +0613 % Genes +0614 for i={'geneMiriams','geneShortNames','proteins'} % Empty strings +0615 model = emptyOrFill(model,i{1},{''},'genes'); +0616 end +0617 for i={'geneComps'} % Ones, assume first compartment +0618 model = emptyOrFill(model,i{1},1,'genes'); +0619 end +0620 % Comps +0621 for i={'compNames'} % Empty strings +0622 model = emptyOrFill(model,i{1},{''},'comps'); +0623 end +0624 for i={'compOutside'} % First comp +0625 model = emptyOrFill(model,i{1},model.comps{1},'comps'); +0626 end +0627 % Single fields are kept, even if empty +0628 % for i={'description','name','version','date','annotation'} +0629 % if isempty(model.(i{1})) +0630 % model = rmfield(model,i{1}); +0631 % end +0632 % end +0633 +0634 if isfield(model,'subSystems') +0635 % If all entries are 1x1, then flatten +0636 if all(cellfun(@(x) numel(x) <= 1, model.subSystems)) +0637 model.subSystems = transpose([model.subSystems{:}]); +0638 end +0639 end +0640 +0641 % Make rxnGeneMat fields and map to the existing model.genes field +0642 if isfield(model,'grRules') +0643 [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules); +0644 model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes)); +0645 [~,geneOrder] = ismember(genes,model.genes); +0646 if any(geneOrder == 0) +0647 error(['The grRules includes the following gene(s), that are not in '... +0648 'the list of model genes: ', genes{~geneOrder}]) 0649 end -0650 end -0651 -0652 if verbose -0653 fprintf(' Done!\n'); -0654 end -0655 end -0656 -0657 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty) -0658 if nargin<5 -0659 keepEmpty=false; -0660 end -0661 if isnumeric(emptyEntry) -0662 emptyCells=isempty(model.(field)); -0663 else -0664 emptyCells=cellfun('isempty',model.(field)); -0665 end -0666 if all(emptyCells) && ~keepEmpty -0667 model = rmfield(model, field); -0668 elseif numel(model.(field))<numel(model.(type)) -0669 model.(field)(end+1:numel(model.(type)),1)=emptyEntry; -0670 end -0671 end -0672 -0673 function model = readFieldValue(model, fieldName, value, pos) -0674 if numel(model.(fieldName))<pos-1 -0675 model.(fieldName)(end+1:pos,1) = {''}; -0676 end -0677 model.(fieldName)(pos,1) = {value}; -0678 end -0679 -0680 function [miriams,miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber) -0681 if isempty(key) -0682 key=miriamKey; -0683 else -0684 miriamKey=key; -0685 end -0686 if ~isempty(value) -0687 miriams(entryNumber,1:3) = {pos, key, strip(value)}; -0688 entryNumber = entryNumber + 1; +0650 model.rxnGeneMat(:,geneOrder) = rxnGeneMat; +0651 end +0652 +0653 % Finalize GECKO model +0654 if isGECKO +0655 % Fill in empty fields and empty entries +0656 for i={'kcat','source','notes','eccodes'} % Even keep empty +0657 model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true); +0658 end +0659 for i={'enzymes','mw','sequence'} +0660 model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true); +0661 end +0662 model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true); +0663 model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true); +0664 % Change string to double +0665 for i={'kcat','mw','concs'} +0666 if isfield(model.ec,i{1}) +0667 model.ec.(i{1}) = str2double(model.ec.(i{1})); +0668 end +0669 end +0670 % Fill rxnEnzMat +0671 rxnIdx = cellfun('isempty', enzStoich(:,1)); +0672 enzStoich(rxnIdx,:) = ''; +0673 rxnIdx = cell2mat(enzStoich(:,1)); +0674 [~,enzIdx] = ismember(enzStoich(:,2),model.ec.enzymes); +0675 coeffs = cell2mat(enzStoich(:,3)); +0676 model.ec.rxnEnzMat = zeros(numel(model.ec.rxns), numel(model.ec.genes)); +0677 linearIndices = sub2ind([numel(model.ec.rxns), numel(model.ec.genes)], rxnIdx, enzIdx); +0678 model.ec.rxnEnzMat(linearIndices) = coeffs; +0679 %Parse ec-codes +0680 if ~isempty(ecGecko) +0681 locs = cell2mat(ecGecko(:,1)); +0682 for i=unique(locs)' +0683 ecGeckoCat=strjoin(ecGecko(locs==i,2),';'); +0684 model.ec.eccodes{i,1}=ecGeckoCat; +0685 end +0686 emptyEc=cellfun('isempty',model.ec.eccodes); +0687 model.ec.eccodes(emptyEc)={''}; +0688 end 0689 end -0690 end +0690 +0691 if verbose +0692 fprintf(' Done!\n'); +0693 end +0694 end +0695 +0696 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty) +0697 if nargin<5 +0698 keepEmpty=false; +0699 end +0700 if isnumeric(emptyEntry) +0701 emptyCells=isempty(model.(field)); +0702 else +0703 emptyCells=cellfun('isempty',model.(field)); +0704 end +0705 if all(emptyCells) && ~keepEmpty +0706 model = rmfield(model, field); +0707 elseif numel(model.(field))<numel(model.(type)) +0708 model.(field)(end+1:numel(model.(type)),1)=emptyEntry; +0709 end +0710 end +0711 +0712 function model = readFieldValue(model, fieldName, value, pos) +0713 if numel(model.(fieldName))<pos-1 +0714 model.(fieldName)(end+1:pos,1) = {''}; +0715 end +0716 model.(fieldName)(pos,1) = {value}; +0717 end +0718 +0719 function [miriams,miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber) +0720 if isempty(key) +0721 key=miriamKey; +0722 else +0723 miriamKey=key; +0724 end +0725 if ~isempty(value) +0726 miriams(entryNumber,1:3) = {pos, key, strip(value)}; +0727 entryNumber = entryNumber + 1; +0728 end +0729 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/solver/optimizeProb.html b/doc/solver/optimizeProb.html index 5fd16c557..328b92faa 100644 --- a/doc/solver/optimizeProb.html +++ b/doc/solver/optimizeProb.html @@ -144,194 +144,203 @@

    SOURCE CODE ^if ~isempty(getCurrentTask) % If run in parallel, then one thread per gurobi -0094 solverparams.Threads=1; -0095 end -0096 solverparams = structUpdate(solverparams,params); -0097 -0098 % Restructering problem according to gurobi format -0099 if isfield(prob, 'csense') -0100 prob.sense = renameparams(prob.csense, {'L','G','E'}, {'<','>','='}); -0101 prob = rmfield(prob, {'csense'}); -0102 end -0103 if isfield(prob, 'osense') -0104 osense = prob.osense; -0105 prob.modelsense = renameparams(num2str(prob.osense), {'1','-1'}, {'min','max'}); -0106 prob = rmfield(prob, {'osense'}); -0107 end -0108 [prob.obj, prob.rhs, prob.vtype] = deal(prob.c, prob.b, prob.vartype); -0109 prob = rmfield(prob, {'c','b','vartype'}); -0110 -0111 resG = gurobi(prob,solverparams); -0112 -0113 try -0114 % Name output fields the same as COBRA does -0115 res.full = resG.x; -0116 res.obj = resG.objval; -0117 res.origStat = resG.status; -0118 if isfield(resG,{'pi','rc'}) -0119 res.dual = -resG.pi*osense; -0120 res.rcost = -resG.rc*osense; -0121 end -0122 if milp && strcmp(resG.status, 'TIME_LIMIT') -0123 % If res has the objval field, it succeeded, regardless of -0124 % time_limit status -0125 resG.status = 'OPTIMAL'; -0126 end -0127 switch resG.status -0128 case 'OPTIMAL' -0129 res.stat = 1; -0130 case 'UNBOUNDED' -0131 res.stat = 2; -0132 otherwise -0133 res.stat = 0; -0134 end -0135 if ~milp -0136 res.vbasis = resG.vbasis; -0137 res.cbasis = resG.cbasis; -0138 else -0139 res.mipgap = resG.mipgap; -0140 end -0141 catch -0142 res.stat = 0; -0143 res.origStat = resG.status; % useful information to have -0144 end -0145 %% Use GLPK using RAVEN-provided binary -0146 case 'glpk' -0147 solverparams.scale = 1; % Auto scaling -0148 %solverparams.tmlim = defaultparams.timeLimit; -0149 solverparams.tolbnd = defaultparams.feasTol; -0150 solverparams.toldj = defaultparams.optTol; -0151 solverparams.tolint = defaultparams.intTol; -0152 solverparams.tolobj = defaultparams.objTol; -0153 solverparams.msglev = 0; % Level of verbosity -0154 solverparams = structUpdate(solverparams,params); -0155 -0156 prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'}); -0157 -0158 if milp -0159 solverparams.tmlim = solverparams.tmlim*10; -0160 solverparams.msglev = 1; % Level of verbosity -0161 disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.') -0162 end -0163 -0164 % Ensure that RAVEN glpk binary is used, return to original -0165 % directory afterwards -0166 [ravenDir,currDir]=findRAVENroot(); -0167 cd(fullfile(ravenDir,'software','GLPKmex')) -0168 [xopt, fmin, errnum, extra] = glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams); -0169 cd(currDir) -0170 -0171 switch errnum % 1 = undefined; 2 = feasible; 3 = infeasible; 4 = no feasible solution; 5 = optimal; 6 = no unbounded solution -0172 case 5 -0173 res.stat = 1; % Optimal -0174 case 2 -0175 res.stat = 2; % Feasible, but not optimal -0176 otherwise -0177 res.stat = 0; -0178 end -0179 res.origStat = errnum; -0180 res.full = xopt; -0181 res.obj = fmin; -0182 res.dual = -extra.lambda*prob.osense; -0183 res.rcost = -extra.redcosts*prob.osense; -0184 %% Use scip -0185 case {'soplex','scip'} % Old 'soplex' option also allowed -0186 [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype); -0187 -0188 % [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts) -0189 % -0190 % Input arguments*: -0191 % H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU]) -0192 % f - linear objective vector -0193 % A - linear constraint matrix (sparse) -0194 % rl - linear constraint lhs -0195 % ru - linear constraint rhs -0196 % lb - decision variable lower bounds -0197 % ub - decision variable upper bounds -0198 % xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary) -0199 % sos - SOS structure with fields type, index and weight (see below) -0200 % qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below) -0201 % nl - Nonlinear Objective and Constraints structure (see below) -0202 % x0 - primal solution -0203 % opts - solver options (see below) -0204 % -0205 % Return arguments: -0206 % x - solution vector -0207 % fval - objective value at the solution -0208 % exitflag - exit status (see below) -0209 % stats - statistics structure -0210 % -0211 % Option Fields (all optional, see also optiset for a list): -0212 % solverOpts - specific SCIP options (list of pairs of parameter names and values) -0213 % maxiter - maximum LP solver iterations -0214 % maxnodes - maximum nodes to explore -0215 % maxtime - maximum execution time [s] -0216 % tolrfun - primal feasibility tolerance -0217 % display - solver display level [0-5] -0218 % probfile - write problem to given file -0219 % presolvedfile - write presolved problem to file -0220 % -0221 % Return Status: -0222 % 0 - Unknown -0223 % 1 - User Interrupted -0224 % 2 - Node Limit Reached -0225 % 3 - Total Node Limit Reached -0226 % 4 - Stall Node Limit Reached -0227 % 5 - Time Limit Reached -0228 % 6 - Memory Limit Reached -0229 % 7 - Gap Limit Reached -0230 % 8 - Solution Limit Reached -0231 % 9 - Solution Improvement Limit Reached -0232 % 10 - Restart Limit Reached -0233 % 11 - Problem Solved to Optimality -0234 % 12 - Problem is Infeasible -0235 % 13 - Problem is Unbounded -0236 % 14 - Problem is Either Infeasible or Unbounded -0237 -0238 res.origStat = exitflag; -0239 res.full = xopt; -0240 res.obj = fval; -0241 -0242 switch exitflag -0243 case 11 -0244 res.stat = 1; -0245 case {1, 5, 6, 7, 8, 9, 10, 13} % If user-interrupted, not sure if optimal -0246 res.stat = 2; -0247 otherwise -0248 res.stat = 0; -0249 end -0250 otherwise -0251 error('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').'); -0252 end -0253 if res.stat>0 -0254 res.full=res.full(1:size(prob.a,2)); -0255 end -0256 end -0257 -0258 function s_merged=structUpdate(s_old,s_new) -0259 %Remove overlapping fields from first struct; -0260 %Obtain all unique names of remaining fields; -0261 %Merge both structs -0262 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new))); -0263 names = [fieldnames(s_merged); fieldnames(s_new)]; -0264 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1); +0093 % If run in parallel, then one thread per gurobi. getCurrentTask +0094 % requires the Parallel Computing Toolbox. +0095 inParallel = false; +0096 if license('test','Distrib_Computing_Toolbox') +0097 try +0098 inParallel = ~isempty(getCurrentTask); +0099 catch +0100 end +0101 end +0102 if inParallel +0103 solverparams.Threads=1; +0104 end +0105 solverparams = structUpdate(solverparams,params); +0106 +0107 % Restructering problem according to gurobi format +0108 if isfield(prob, 'csense') +0109 prob.sense = renameparams(prob.csense, {'L','G','E'}, {'<','>','='}); +0110 prob = rmfield(prob, {'csense'}); +0111 end +0112 if isfield(prob, 'osense') +0113 osense = prob.osense; +0114 prob.modelsense = renameparams(num2str(prob.osense), {'1','-1'}, {'min','max'}); +0115 prob = rmfield(prob, {'osense'}); +0116 end +0117 [prob.obj, prob.rhs, prob.vtype] = deal(prob.c, prob.b, prob.vartype); +0118 prob = rmfield(prob, {'c','b','vartype'}); +0119 +0120 resG = gurobi(prob,solverparams); +0121 +0122 try +0123 % Name output fields the same as COBRA does +0124 res.full = resG.x; +0125 res.obj = resG.objval; +0126 res.origStat = resG.status; +0127 if isfield(resG,{'pi','rc'}) +0128 res.dual = -resG.pi*osense; +0129 res.rcost = -resG.rc*osense; +0130 end +0131 if milp && strcmp(resG.status, 'TIME_LIMIT') +0132 % If res has the objval field, it succeeded, regardless of +0133 % time_limit status +0134 resG.status = 'OPTIMAL'; +0135 end +0136 switch resG.status +0137 case 'OPTIMAL' +0138 res.stat = 1; +0139 case 'UNBOUNDED' +0140 res.stat = 2; +0141 otherwise +0142 res.stat = 0; +0143 end +0144 if ~milp +0145 res.vbasis = resG.vbasis; +0146 res.cbasis = resG.cbasis; +0147 else +0148 res.mipgap = resG.mipgap; +0149 end +0150 catch +0151 res.stat = 0; +0152 res.origStat = resG.status; % useful information to have +0153 end +0154 %% Use GLPK using RAVEN-provided binary +0155 case 'glpk' +0156 solverparams.scale = 1; % Auto scaling +0157 %solverparams.tmlim = defaultparams.timeLimit; +0158 solverparams.tolbnd = defaultparams.feasTol; +0159 solverparams.toldj = defaultparams.optTol; +0160 solverparams.tolint = defaultparams.intTol; +0161 solverparams.tolobj = defaultparams.objTol; +0162 solverparams.msglev = 0; % Level of verbosity +0163 solverparams = structUpdate(solverparams,params); +0164 +0165 prob.csense = renameparams(prob.csense, {'L','G','E'}, {'U','L','S'}); +0166 +0167 if milp +0168 solverparams.tmlim = solverparams.tmlim*10; +0169 solverparams.msglev = 1; % Level of verbosity +0170 disp('Issues have been observed when using GLPK for MILP solving. Be advised to carefully observe the results, or us another solver.') +0171 end +0172 +0173 % Ensure that RAVEN glpk binary is used, return to original +0174 % directory afterwards +0175 [ravenDir,currDir]=findRAVENroot(); +0176 cd(fullfile(ravenDir,'software','GLPKmex')) +0177 [xopt, fmin, errnum, extra] = glpk(prob.c, prob.A, prob.b, prob.lb, prob.ub, prob.csense, prob.vartype, prob.osense, solverparams); +0178 cd(currDir) +0179 +0180 switch errnum % 1 = undefined; 2 = feasible; 3 = infeasible; 4 = no feasible solution; 5 = optimal; 6 = no unbounded solution +0181 case 5 +0182 res.stat = 1; % Optimal +0183 case 2 +0184 res.stat = 2; % Feasible, but not optimal +0185 otherwise +0186 res.stat = 0; +0187 end +0188 res.origStat = errnum; +0189 res.full = xopt; +0190 res.obj = fmin; +0191 res.dual = -extra.lambda*prob.osense; +0192 res.rcost = -extra.redcosts*prob.osense; +0193 %% Use scip +0194 case {'soplex','scip'} % Old 'soplex' option also allowed +0195 [xopt,fval,exitflag] = scip([], prob.c, prob.A,-prob.b, prob.b, prob.lb, prob.ub, prob.vartype); +0196 +0197 % [x,fval,exitflag,stats] = scip(H, f, A, rl, ru, lb, ub, xtype, sos, qc, nl, x0, opts) +0198 % +0199 % Input arguments*: +0200 % H - quadratic objective matrix (sparse, optional [NOT TRIL / TRIU]) +0201 % f - linear objective vector +0202 % A - linear constraint matrix (sparse) +0203 % rl - linear constraint lhs +0204 % ru - linear constraint rhs +0205 % lb - decision variable lower bounds +0206 % ub - decision variable upper bounds +0207 % xtype - string of variable integrality ('c' continuous, 'i' integer, 'b' binary) +0208 % sos - SOS structure with fields type, index and weight (see below) +0209 % qc - Quadratic Constraints structure with fields Q, l, qrl and qru (see below) +0210 % nl - Nonlinear Objective and Constraints structure (see below) +0211 % x0 - primal solution +0212 % opts - solver options (see below) +0213 % +0214 % Return arguments: +0215 % x - solution vector +0216 % fval - objective value at the solution +0217 % exitflag - exit status (see below) +0218 % stats - statistics structure +0219 % +0220 % Option Fields (all optional, see also optiset for a list): +0221 % solverOpts - specific SCIP options (list of pairs of parameter names and values) +0222 % maxiter - maximum LP solver iterations +0223 % maxnodes - maximum nodes to explore +0224 % maxtime - maximum execution time [s] +0225 % tolrfun - primal feasibility tolerance +0226 % display - solver display level [0-5] +0227 % probfile - write problem to given file +0228 % presolvedfile - write presolved problem to file +0229 % +0230 % Return Status: +0231 % 0 - Unknown +0232 % 1 - User Interrupted +0233 % 2 - Node Limit Reached +0234 % 3 - Total Node Limit Reached +0235 % 4 - Stall Node Limit Reached +0236 % 5 - Time Limit Reached +0237 % 6 - Memory Limit Reached +0238 % 7 - Gap Limit Reached +0239 % 8 - Solution Limit Reached +0240 % 9 - Solution Improvement Limit Reached +0241 % 10 - Restart Limit Reached +0242 % 11 - Problem Solved to Optimality +0243 % 12 - Problem is Infeasible +0244 % 13 - Problem is Unbounded +0245 % 14 - Problem is Either Infeasible or Unbounded +0246 +0247 res.origStat = exitflag; +0248 res.full = xopt; +0249 res.obj = fval; +0250 +0251 switch exitflag +0252 case 11 +0253 res.stat = 1; +0254 case {1, 5, 6, 7, 8, 9, 10, 13} % If user-interrupted, not sure if optimal +0255 res.stat = 2; +0256 otherwise +0257 res.stat = 0; +0258 end +0259 otherwise +0260 error('RAVEN solver not defined or unknown. Try using setRavenSolver(''solver'').'); +0261 end +0262 if res.stat>0 +0263 res.full=res.full(1:size(prob.a,2)); +0264 end 0265 end 0266 -0267 function paramlist = renameparams(paramlist,old,new) -0268 if ~iscell(paramlist) -0269 wasNoCell = true; -0270 paramlist={paramlist}; -0271 else -0272 wasNoCell = false; -0273 end -0274 for i=1:numel(old) -0275 paramlist = regexprep(paramlist,old{i},new{i}); -0276 end -0277 if wasNoCell -0278 paramlist=paramlist{1}; -0279 end -0280 end +0267 function s_merged=structUpdate(s_old,s_new) +0268 %Remove overlapping fields from first struct; +0269 %Obtain all unique names of remaining fields; +0270 %Merge both structs +0271 s_merged = rmfield(s_old, intersect(fieldnames(s_old), fieldnames(s_new))); +0272 names = [fieldnames(s_merged); fieldnames(s_new)]; +0273 s_merged = cell2struct([struct2cell(s_merged); struct2cell(s_new)], names, 1); +0274 end +0275 +0276 function paramlist = renameparams(paramlist,old,new) +0277 if ~iscell(paramlist) +0278 wasNoCell = true; +0279 paramlist={paramlist}; +0280 else +0281 wasNoCell = false; +0282 end +0283 for i=1:numel(old) +0284 paramlist = regexprep(paramlist,old{i},new{i}); +0285 end +0286 if wasNoCell +0287 paramlist=paramlist{1}; +0288 end +0289 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/testing/unit_tests/modelCurationTests.html b/doc/testing/unit_tests/modelCurationTests.html index dde7fe80c..0846789f0 100644 --- a/doc/testing/unit_tests/modelCurationTests.html +++ b/doc/testing/unit_tests/modelCurationTests.html @@ -42,7 +42,7 @@

    CROSS-REFERENCE INFORMATION ^ <h2><a name=SUBFUNCTIONS ^

    +
  • function addExchangeRxnsTest(testCase)
  • function addTransportTest(testCase)
  • function addGenesRavenTest(testCase)
  • function removeGenesTest(testCase)
  • function addMetsTest(testCase)
  • function addMets_oneCompTest(testCase)
  • function removeMetsTest(testCase)
  • function addRxnsTest(testCase)
  • function removeReactionsTest(testCase)
  • function getMetsInCompTest(testCase)
  • function getRxnsInCompTest(testCase)
  • function getAllRxnsFromGenesTest(testCase)
  • function changeRxnsTest(testCase)
  • function removeBadRxnsTest(testCase)
  • function changeGrRulesTest(testCase)
  • function copyToCompsTest(testCase)
  • function setExchangeBoundsTest(testCase)
  • function addRxnsGenesMetsTest(testCase)
  • function replaceMetsOutputsTest(testCase)
  • function replaceMetsSharedReactionTest(testCase)
  • function replaceMetsNamesTest(testCase)
  • function replaceMetsNamesRepeatedReplacementTest(testCase)

    • SOURCE CODE ^

    0001 %run this test case with the command
@@ -537,7 +537,95 @@ 
    SOURCE CODE ^end
    +0493 end +0494 +0495 function replaceMetsOutputsTest(testCase) +0496 sourceDir = fileparts(which(mfilename)); +0497 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0498 +0499 %Add a duplicate of atp_c and a copy of ATPM that uses it +0500 metsToAdd.mets={'atpDup_c'}; +0501 metsToAdd.metNames={'atp duplicate'}; +0502 metsToAdd.compartments={'c'}; +0503 evalc('model=addMets(model,metsToAdd);'); +0504 rxnsToAdd.rxns={'ATPM_dup'}; +0505 rxnsToAdd.equations={'atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +0506 evalc('model=addRxns(model,rxnsToAdd,1);'); +0507 +0508 [modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'atpDup_c','atp_c',false,true); +0509 +0510 verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +0511 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +0512 verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +0513 verifyEqual(testCase,idxDuplRxns,numel(model.rxns)) +0514 end +0515 +0516 function replaceMetsSharedReactionTest(testCase) +0517 %The replaced metabolite occurs in a reaction that already contains the +0518 %replacement metabolite; the stoichiometries are summed. +0519 sourceDir = fileparts(which(mfilename)); +0520 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0521 +0522 metsToAdd.mets={'atpDup_c'}; +0523 metsToAdd.metNames={'atp duplicate'}; +0524 metsToAdd.compartments={'c'}; +0525 evalc('model=addMets(model,metsToAdd);'); +0526 rxnsToAdd.rxns={'ATPshared'}; +0527 rxnsToAdd.equations={'atp_c + atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +0528 evalc('model=addRxns(model,rxnsToAdd,1);'); +0529 +0530 modelNew=replaceMets(model,'atpDup_c','atp_c',false,true); +0531 +0532 rxnIdx=strcmp(modelNew.rxns,'ATPshared'); +0533 verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),rxnIdx)),-2) +0534 verifyEqual(testCase,numel(modelNew.mets),numel(model.mets)-1) +0535 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +0536 end +0537 +0538 function replaceMetsNamesTest(testCase) +0539 %The replacement metabolite has a higher index than the to-be-replaced +0540 %metabolite, and the removed duplicate reaction is not the last reaction. +0541 sourceDir = fileparts(which(mfilename)); +0542 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0543 +0544 metsToAdd.mets={'atpNew_c'}; +0545 metsToAdd.metNames={'ATPnew'}; +0546 metsToAdd.compartments={'c'}; +0547 evalc('model=addMets(model,metsToAdd);'); +0548 rxnsToAdd.rxns={'ATPM_dup','ATPM_other'}; +0549 rxnsToAdd.equations={'atpNew_c + h2o_c => adp_c + h_c + pi_c','h2o_c => h_c + pi_c'}; +0550 evalc('model=addRxns(model,rxnsToAdd,1);'); +0551 +0552 [modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'ATP','ATPnew'); +0553 +0554 verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +0555 verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +0556 verifyEqual(testCase,idxDuplRxns,numel(model.rxns)-1) +0557 verifyEqual(testCase,model.rxns(idxDuplRxns),removedRxns) +0558 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp_c'))) +0559 end +0560 +0561 function replaceMetsNamesRepeatedReplacementTest(testCase) +0562 %The replacement name occurs twice in the same compartment; the +0563 %stoichiometry is added once. +0564 sourceDir = fileparts(which(mfilename)); +0565 load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); +0566 +0567 metsToAdd.mets={'atp2_c','atpNew_c'}; +0568 metsToAdd.metNames={'ATP2','ATPnew'}; +0569 metsToAdd.compartments={'c','c'}; +0570 evalc('model=addMets(model,metsToAdd);'); +0571 model.metNames{strcmp(model.mets,'atp2_c')}='ATP'; +0572 rxnsToAdd.rxns={'rNew'}; +0573 rxnsToAdd.equations={'atpNew_c => adp_c'}; +0574 evalc('model=addRxns(model,rxnsToAdd,1);'); +0575 +0576 modelNew=replaceMets(model,'ATPnew','ATP'); +0577 +0578 verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),strcmp(modelNew.rxns,'rNew'))),-1) +0579 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpNew_c'))) +0580 verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp2_c'))) +0581 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/installation/checkInstallation.m b/installation/checkInstallation.m index de8c12363..40ef456fe 100755 --- a/installation/checkInstallation.m +++ b/installation/checkInstallation.m @@ -309,8 +309,37 @@ [~,res]=evalc("runtests('hmmerTests.m')"); res=interpretResults(res); if res==false - fprintf([' This is essential to run getKEGGModelFromHomology()\n'... - ' when using a FASTA file as input\n']) + fprintf([' HMMER is only required by getKEGGModelForOrganism() when\n'... + ' using a FASTA file as input. If you do not use this function,\n'... + ' this failure can be ignored.\n']) + % On Windows, HMMER runs through Windows Subsystem for Linux + if ispc + [status,~]=system('where wsl'); + if status~=0 + printOrange([' MATLAB cannot find wsl.exe. On Windows, HMMER runs through WSL.\n'... + ' Install it with "wsl --install" in a Command Prompt, then restart\n'... + ' MATLAB (it inherits the system PATH at startup) and rerun checkInstallation\n']) + else + [status,~]=system('wsl -e true'); + if status~=0 + printOrange([' MATLAB finds wsl.exe, but cannot start a Linux distribution\n'... + ' through it. Run "wsl -l -v" in a Command Prompt to check that a\n'... + ' distribution is installed and set up, then restart MATLAB and\n'... + ' rerun checkInstallation\n']) + else + try + [status,~]=system(['wsl "' getWSLpath(fullfile(ravenDir,'software','hmmer','hmmsearch')) '" -h']); + catch + status=1; + end + if status~=0 + printOrange([' WSL works, but the RAVEN HMMER binary cannot be run in it.\n'... + ' Check that the files in software/hmmer are present and that\n'... + ' the drive holding RAVEN is accessible from WSL\n']) + end + end + end + end end if developMode diff --git a/io/importModel.m b/io/importModel.m index dae96063c..b8fdade73 100755 --- a/io/importModel.m +++ b/io/importModel.m @@ -739,11 +739,22 @@ end end if isfield(modelSBML,'notes') - startString=strfind(modelSBML.notes,'xhtml">'); - endString=strfind(modelSBML.notes,''); - if any(startString) && any(endString) - model.annotation.note=modelSBML.notes(startString(1)+7:endString-1); - model.annotation.note=regexprep(model.annotation.note,'|',''); + %Two shapes reach here. RAVEN's own exportModel writes the note as bare + %text inside a ; a writer that serialises a notes dictionary emits + %one '

    key: value

    ' per entry and no at all, carrying the note + %under its own key. Take the labelled entry when there is one, so the + %other entries (version, and whatever else the dictionary held) stay out + %of the note, and fall back to the text otherwise. + note=parseNote(modelSBML.notes,'note'); + if isempty(note) + startString=strfind(modelSBML.notes,'xhtml">'); + endString=strfind(modelSBML.notes,''); + if any(startString) && any(endString) + note=modelSBML.notes(startString(1)+7:endString-1); + end + end + if ~isempty(note) + model.annotation.note=regexprep(note,'|',''); model.annotation.note=strtrim(model.annotation.note); if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML') model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel diff --git a/io/readYAMLmodel.m b/io/readYAMLmodel.m index 1eb65f3b6..ea353f39f 100755 --- a/io/readYAMLmodel.m +++ b/io/readYAMLmodel.m @@ -148,7 +148,7 @@ tline_value = line_value{i}; % import different sections switch tline_raw - case '- metaData:' + case {'- metaData:', '- metaData: !!omap'} section = 1; if verbose fprintf('\t%d\n', section); @@ -168,7 +168,7 @@ end pos=0; continue - case '- genes:' + case {'- genes:', '- genes: []'} section = 4; if verbose fprintf('\t%d\n', section); @@ -234,7 +234,7 @@ model.annotation.email = tline_value; case 'organization' model.annotation.organization = tline_value; - case 'geckoLight' + case {'geckoLight', 'gecko_light'} if strcmp(tline_value,'true') model.ec.geckoLight = true; end @@ -267,8 +267,12 @@ model = readFieldValue(model, 'inchis', tline_value, pos); readList=''; miriamKey=''; case 'smiles' - model = readFieldValue(model, 'metSmiles', tline_value, pos); - readList=''; miriamKey=''; + if isempty(tline_value) + readList = 'smilesInAnnotation'; + else + model = readFieldValue(model, 'metSmiles', tline_value, pos); + readList=''; miriamKey=''; + end case 'deltaG' model = readFieldValue(model, 'metDeltaG', tline_value, pos); readList=''; miriamKey=''; @@ -279,6 +283,16 @@ readList = 'annotation'; otherwise switch readList + case 'smilesInAnnotation' + %As for ec-code on a reaction: a bare list item + %continues the list, anything else ends it. + if isempty(tline_key) + model = readFieldValue(model, 'metSmiles', ... + regexprep(tline_value,'^ +- "?(.*)"?$','$1'), pos); + else + readList = 'annotation'; + [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); + end case 'annotation' [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo); otherwise @@ -309,7 +323,7 @@ case 'gene_reaction_rule' model = readFieldValue(model, 'grRules', tline_value, pos); readList=''; miriamKey=''; - case 'rxnNotes' + case {'rxnNotes', 'notes'} model = readFieldValue(model, 'rxnNotes', tline_value, pos); readList=''; miriamKey=''; case 'rxnFrom' @@ -343,6 +357,19 @@ end case 'metabolites' readList = 'equation'; + case 'ec-code' + %EC numbers emitted inside the annotation block rather than + %under the legacy top-level 'eccodes' key, which is still + %accepted above. A block list follows when the value is + %empty; it is read under its own state so that gathering can + %resume for the rest of the annotation block once the list + %ends, instead of swallowing every later entry as an EC code. + if isempty(tline_value) + readList = 'eccodesInAnnotation'; + else + eccodes(ecCodeNo,1:2)={pos,tline_value}; + ecCodeNo=ecCodeNo+1; + end case 'annotation' readList = 'annotation'; @@ -351,6 +378,18 @@ case 'eccodes' eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; ecCodeNo=ecCodeNo+1; + case 'eccodesInAnnotation' + %A bare list item ("- 1.1.2.4") has no colon, so + %tline_key is empty and the list continues. Anything + %else is the annotation block's next entry, so the + %list has ended. + if isempty(tline_key) + eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; + ecCodeNo=ecCodeNo+1; + else + readList = 'annotation'; + [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo); + end case 'subsystem' subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')}; subSysNo=subSysNo+1; diff --git a/readme/RAVEN_structure_fields.csv b/readme/RAVEN_structure_fields.csv index f471abf7b..4c95e19cb 100644 --- a/readme/RAVEN_structure_fields.csv +++ b/readme/RAVEN_structure_fields.csv @@ -49,7 +49,7 @@ geneMiriams;;Struct;Structure with MIRIAM information about the genes;Non-mandat ;value;Column Cell Array of Strings;The ID for particular MIRIAM entry;Non-mandatory field;P12345; geneShortNames;;Column Cell Array of Strings;Gene alternative names (e.g. ERG10);Non-mandatory field;ERG10; metNames;;Column Cell Array of Strings;Metabolite names;Non-mandatory field;alpha-D-glucose; -metComps;;Column Vector of Doubles;Compartments for metabolites;Mandatory field;1;Shows the index of affiliated compartment in comps, not the compartment abbreviation +metComps;;Column Vector of Integers;Compartments for metabolites;Mandatory field;1;Written as an integer so a reader that converts the 1-based index into a position gets one; MATLAB stores an integer-valued double narrowed to an integer type anyway, and both classes index comps identically. inchis;;Column Cell Array of Strings;InChI-codes for metabolites;Non-mandatory field;1/C6H12O6/c7-1-2-3(8)4(9)5(10)6(11)12-2/h2-11H,1H2/t2-,3-,4+,5-,6+/m1/s1; metFormulas;;Column Cell Array of Strings;Chemical formulas for metabolites;Non-mandatory field;C6H12O6; metMiriams;;Struct;Structure with MIRIAM information about the metabolites;Non-mandatory field;;If there are multiple miriams, 'name' indexes must match 'value' indexes diff --git a/solver/optimizeProb.m b/solver/optimizeProb.m index 3430c45d3..bcb137e98 100755 --- a/solver/optimizeProb.m +++ b/solver/optimizeProb.m @@ -90,7 +90,16 @@ solverparams.FeasibilityTol = defaultparams.feasTol; solverparams.OptimalityTol = defaultparams.optTol; solverparams.Presolve = 2; - if ~isempty(getCurrentTask) % If run in parallel, then one thread per gurobi + % If run in parallel, then one thread per gurobi. getCurrentTask + % requires the Parallel Computing Toolbox. + inParallel = false; + if license('test','Distrib_Computing_Toolbox') + try + inParallel = ~isempty(getCurrentTask); + catch + end + end + if inParallel solverparams.Threads=1; end solverparams = structUpdate(solverparams,params); diff --git a/testing/unit_tests/modelCurationTests.m b/testing/unit_tests/modelCurationTests.m index 62e3dc962..bd6f80900 100755 --- a/testing/unit_tests/modelCurationTests.m +++ b/testing/unit_tests/modelCurationTests.m @@ -491,3 +491,91 @@ function addRxnsGenesMetsTest(testCase) verifyEqual(testCase,modelNew,modelManual) end + +function replaceMetsOutputsTest(testCase) +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +%Add a duplicate of atp_c and a copy of ATPM that uses it +metsToAdd.mets={'atpDup_c'}; +metsToAdd.metNames={'atp duplicate'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPM_dup'}; +rxnsToAdd.equations={'atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +[modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'atpDup_c','atp_c',false,true); + +verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +verifyEqual(testCase,idxDuplRxns,numel(model.rxns)) +end + +function replaceMetsSharedReactionTest(testCase) +%The replaced metabolite occurs in a reaction that already contains the +%replacement metabolite; the stoichiometries are summed. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atpDup_c'}; +metsToAdd.metNames={'atp duplicate'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPshared'}; +rxnsToAdd.equations={'atp_c + atpDup_c + h2o_c => adp_c + h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +modelNew=replaceMets(model,'atpDup_c','atp_c',false,true); + +rxnIdx=strcmp(modelNew.rxns,'ATPshared'); +verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),rxnIdx)),-2) +verifyEqual(testCase,numel(modelNew.mets),numel(model.mets)-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpDup_c'))) +end + +function replaceMetsNamesTest(testCase) +%The replacement metabolite has a higher index than the to-be-replaced +%metabolite, and the removed duplicate reaction is not the last reaction. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atpNew_c'}; +metsToAdd.metNames={'ATPnew'}; +metsToAdd.compartments={'c'}; +evalc('model=addMets(model,metsToAdd);'); +rxnsToAdd.rxns={'ATPM_dup','ATPM_other'}; +rxnsToAdd.equations={'atpNew_c + h2o_c => adp_c + h_c + pi_c','h2o_c => h_c + pi_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +[modelNew,removedRxns,idxDuplRxns]=replaceMets(model,'ATP','ATPnew'); + +verifyEqual(testCase,numel(modelNew.rxns),numel(model.rxns)-1) +verifyEqual(testCase,removedRxns,{'ATPM_dup'}) +verifyEqual(testCase,idxDuplRxns,numel(model.rxns)-1) +verifyEqual(testCase,model.rxns(idxDuplRxns),removedRxns) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp_c'))) +end + +function replaceMetsNamesRepeatedReplacementTest(testCase) +%The replacement name occurs twice in the same compartment; the +%stoichiometry is added once. +sourceDir = fileparts(which(mfilename)); +load(fullfile(sourceDir,'test_data','ecoli_textbook.mat'), 'model'); + +metsToAdd.mets={'atp2_c','atpNew_c'}; +metsToAdd.metNames={'ATP2','ATPnew'}; +metsToAdd.compartments={'c','c'}; +evalc('model=addMets(model,metsToAdd);'); +model.metNames{strcmp(model.mets,'atp2_c')}='ATP'; +rxnsToAdd.rxns={'rNew'}; +rxnsToAdd.equations={'atpNew_c => adp_c'}; +evalc('model=addRxns(model,rxnsToAdd,1);'); + +modelNew=replaceMets(model,'ATPnew','ATP'); + +verifyEqual(testCase,full(modelNew.S(strcmp(modelNew.mets,'atp_c'),strcmp(modelNew.rxns,'rNew'))),-1) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atpNew_c'))) +verifyEmpty(testCase,find(strcmp(modelNew.mets,'atp2_c'))) +end