From 162eeaf88d95156ba20f2be831b67c8c717bde42 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 7 Oct 2026 16:13:06 +0200 Subject: [PATCH 1/3] fix: read RAVEN 3 model files, and tolerate a non-double rev RAVEN 3 and raven-toolbox write YAML this version cannot read. The metaData map is tagged (`- metaData: !!omap`), which the exact-match case missed, so the whole block was skipped and the model came back with no id, name, version or date. A reaction's free text moved from `rxnNotes` to `notes`, which aborted the read outright with "Unknown entry in yaml file". EC codes and SMILES moved inside the per-entry `annotation` block, where they were gathered as MIRIAMs instead of filling eccodes and metSmiles. An empty gene list is written inline, and the GECKO flag is spelled `gecko_light`. All of these are accepted on read now; what this version writes is unchanged, so its own output stays valid for older RAVEN. importModel also only found a model note inside a , which RAVEN's exportModel writes but a writer serialising a notes dictionary does not; such a writer emits one "

key: value

" per entry, so the labelled entry is taken when there is one and the text otherwise. mergeLinear combines rev with a sparse logical, which MATLAB refuses for an integer class ("Sparse integer array arithmetic operations are not supported"). A model built outside RAVEN can carry an integer rev, so it is coerced once on entry rather than failing several steps into ftINIT. metComps is documented as an integer column: a reader that turns the 1-based index into a position needs an integer, and MATLAB stores an integer-valued double narrowed to an integer type in any case. --- INIT/mergeLinear.m | 6 ++++ io/importModel.m | 21 ++++++++++--- io/readYAMLmodel.m | 51 +++++++++++++++++++++++++++---- readme/RAVEN_structure_fields.csv | 2 +- 4 files changed, 68 insertions(+), 12 deletions(-) diff --git a/INIT/mergeLinear.m b/INIT/mergeLinear.m index 1378353b2..be642369f 100644 --- a/INIT/mergeLinear.m +++ b/INIT/mergeLinear.m @@ -39,6 +39,12 @@ reducedModel.geneComps=[]; end +%rev is combined with a sparse logical below ('S(...).' > 0 | rev'), which +%MATLAB refuses for an integer class: "Sparse integer array arithmetic +%operations are not supported". A model built outside RAVEN can carry an +%integer rev, so coerce it once here rather than failing several steps in. +reducedModel.rev = double(reducedModel.rev); + nextGroupId = 1; origRxnIds = reducedModel.rxns; groupIds = zeros(numel(reducedModel.rxns),1); 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 From 076b6e213b405403cf03ede00111b04387f30d05 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 7 Oct 2026 17:02:45 +0200 Subject: [PATCH 2/3] docs: regenerate HTML for the changed functions updateDocumentation regenerates the whole doc/ tree, but the committed files were generated on Windows and every other one differs only in the path separator of its index link and title (external\kegg vs external/kegg). Only the three functions this branch changes are updated, so the diff stays about the change rather than about the platform it was regenerated on. --- doc/INIT/mergeLinear.html | 348 +++++++-------- doc/io/importModel.html | 807 +++++++++++++++++----------------- doc/io/readYAMLmodel.html | 879 ++++++++++++++++++++------------------ 3 files changed, 1045 insertions(+), 989 deletions(-) diff --git a/doc/INIT/mergeLinear.html b/doc/INIT/mergeLinear.html index 16467e80e..2c5c663a5 100644 --- a/doc/INIT/mergeLinear.html +++ b/doc/INIT/mergeLinear.html @@ -102,180 +102,186 @@

SOURCE CODE ^end 0041 -0042 nextGroupId = 1; -0043 origRxnIds = reducedModel.rxns; -0044 groupIds = zeros(numel(reducedModel.rxns),1); -0045 reversedRxns = false(numel(reducedModel.rxns),1); -0046 -0047 %Loop through and iteratively group linear reactions -0048 while 1 -0049 %Get the banned reaction indexes. Note that the indexes will change -0050 %in each iteration, but the names will not as they won't be merged -0051 %with any other reaction -0052 bannedIndexes=getIndexes(reducedModel,noMergeRxns,'rxns'); -0053 -0054 %Select all metabolites that are only present as reactants/products -0055 %in one reaction -0056 twoNonZero = find(sum(reducedModel.S ~= 0, 2) == 2); -0057 -0058 mergedSome=false; +0042 %rev is combined with a sparse logical below ('S(...).' > 0 | rev'), which +0043 %MATLAB refuses for an integer class: "Sparse integer array arithmetic +0044 %operations are not supported". A model built outside RAVEN can carry an +0045 %integer rev, so coerce it once here rather than failing several steps in. +0046 reducedModel.rev = double(reducedModel.rev); +0047 +0048 nextGroupId = 1; +0049 origRxnIds = reducedModel.rxns; +0050 groupIds = zeros(numel(reducedModel.rxns),1); +0051 reversedRxns = false(numel(reducedModel.rxns),1); +0052 +0053 %Loop through and iteratively group linear reactions +0054 while 1 +0055 %Get the banned reaction indexes. Note that the indexes will change +0056 %in each iteration, but the names will not as they won't be merged +0057 %with any other reaction +0058 bannedIndexes=getIndexes(reducedModel,noMergeRxns,'rxns'); 0059 -0060 %Loop through each of them and see if the reactions should be -0061 %merged -0062 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))); -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; +0060 %Select all metabolites that are only present as reactants/products +0061 %in one reaction +0062 twoNonZero = find(sum(reducedModel.S ~= 0, 2) == 2); +0063 +0064 mergedSome=false; +0065 +0066 %Loop through each of them and see if the reactions should be +0067 %merged +0068 for i=1:numel(twoNonZero) +0069 involvedRxns=find(reducedModel.S(twoNonZero(i),:)); +0070 %Check that we can have one positive and one negative +0071 pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns)); +0072 neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns)); +0073 +0074 +0075 %Check so that one or both of the reactions haven't been merged +0076 %already +0077 if numel(involvedRxns)==2 && isempty(intersect(bannedIndexes,involvedRxns)) && pos >= 1 && neg >= 1 +0078 %first, take care of a special case: If the first reaction is producing the metabolite and if it is reversible, +0079 %and the second is also producing it and is not reversible, change the order - the code below will not work otherwise +0080 if reducedModel.rev(involvedRxns(1)) && (~reducedModel.rev(involvedRxns(2))) && ... +0081 (reducedModel.S(twoNonZero(i),involvedRxns(1)) > 0) && (reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0) +0082 involvedRxns = flip(involvedRxns); +0083 end +0084 +0085 %first make sure the first reaction is producing the metabolite +0086 if reducedModel.S(twoNonZero(i),involvedRxns(1)) < 0 +0087 %it is not producing the metabolite - fix that +0088 %first choice: use the second reaction as producer if it is producing +0089 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 +0090 involvedRxns = flip(involvedRxns);%make the second the first +0091 else +0092 %now we know that the second reaction is not producing, so we can safely try to make the first a producer +0093 if reducedModel.rev(involvedRxns(1)) == 1 +0094 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(1), groupIds, reversedRxns); +0095 else %ok, finally try to flip the second reaction +0096 if reducedModel.rev(involvedRxns(2)) == 1 +0097 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); +0098 involvedRxns = flip(involvedRxns);%make the second the first +0099 else +0100 error('We should never end up here!'); +0101 end +0102 end +0103 end +0104 end +0105 %Now, make sure the second rxn is a consumer +0106 if reducedModel.S(twoNonZero(i),involvedRxns(2)) > 0 +0107 if reducedModel.rev(involvedRxns(2)) == 1 +0108 [reducedModel,reversedRxns] = flipRxn(reducedModel, involvedRxns(2), groupIds, reversedRxns); +0109 else +0110 error('We should never end up here!'); +0111 end +0112 end +0113 +0114 %Calculate how many times the second reaction has to be +0115 %multiplied before being merged with the first +0116 stoichRatio=abs(reducedModel.S(twoNonZero(i),involvedRxns(1))/reducedModel.S(twoNonZero(i),involvedRxns(2))); 0117 -0118 %This is to prevent numerical issues. It should be 0 -0119 %already -0120 reducedModel.S(twoNonZero(i),involvedRxns(1))=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; -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); +0118 %Add the second to the first +0119 reducedModel.S(:,involvedRxns(1))=reducedModel.S(:,involvedRxns(1))+reducedModel.S(:,involvedRxns(2))*stoichRatio; +0120 +0121 %Clear the second reaction +0122 reducedModel.S(:,involvedRxns(2))=0; +0123 +0124 %This is to prevent numerical issues. It should be 0 +0125 %already +0126 reducedModel.S(twoNonZero(i),involvedRxns(1))=0; +0127 +0128 %At this point the second reaction is certain to be deleted +0129 %in a later step and can therefore be ignored +0130 +0131 %Recalculate the bounds for the new reaction. This can be +0132 %problematic since the scale of the bounds may change +0133 %dramatically. Let the most constraining reaction determine +0134 %the new bound +0135 lb1=reducedModel.lb(involvedRxns(1)); +0136 lb2=reducedModel.lb(involvedRxns(2)); +0137 ub1=reducedModel.ub(involvedRxns(1)); +0138 ub2=reducedModel.ub(involvedRxns(2)); +0139 +0140 if lb2~=-inf +0141 reducedModel.lb(involvedRxns(1))=max(lb1,lb2/stoichRatio); +0142 end +0143 if ub2~=inf +0144 reducedModel.ub(involvedRxns(1))=min(ub1,ub2/stoichRatio); +0145 end +0146 +0147 %take care of the .rev flag - it could be that the combined rxn changes from rev to irrev +0148 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 +0149 +0150 %Then recalculate the objective coefficient. The resulting +0151 %coefficient is the weighted sum of the previous +0152 reducedModel.c(involvedRxns(1))=reducedModel.c(involvedRxns(1))+reducedModel.c(involvedRxns(2))*stoichRatio; +0153 +0154 %store which reactions that have been merged +0155 rxnInd1 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(1)))); +0156 rxnInd2 = find(strcmp(origRxnIds, reducedModel.rxns(involvedRxns(2)))); +0157 grpId = max(groupIds(rxnInd1),groupIds(rxnInd2)); +0158 if grpId == 0 +0159 grpId = nextGroupId; +0160 nextGroupId = nextGroupId + 1; +0161 end +0162 +0163 if groupIds(rxnInd1) ~= grpId +0164 if groupIds(rxnInd1) == 0 +0165 %not merged before, just set the group id +0166 groupIds(rxnInd1) = grpId; +0167 else +0168 %merged before - all rxns with the same group id should be changed +0169 groupIds(groupIds == groupIds(rxnInd1)) = grpId; +0170 end +0171 end +0172 if groupIds(rxnInd2) ~= grpId +0173 if groupIds(rxnInd2) == 0 +0174 %not merged before, just set the group id +0175 groupIds(rxnInd2) = grpId; +0176 else +0177 %merged before - all rxns with the same group id should be changed +0178 groupIds(groupIds == groupIds(rxnInd2)) = grpId; +0179 end +0180 end +0181 +0182 %Iterate again +0183 mergedSome=true; +0184 end +0185 end +0186 +0187 %All possible reactions merged +0188 if mergedSome==false +0189 break; +0190 end 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 +0192 %Now delete all reactions that involve no metabolites +0193 I=find(sum(reducedModel.S~=0,1)==0); +0194 +0195 %Remove reactions +0196 reducedModel=removeReactions(reducedModel,I); +0197 +0198 %Remove metabolites +0199 notInUse=sum(reducedModel.S~=0,2)==0; +0200 reducedModel=removeMets(reducedModel,notInUse); +0201 end +0202 +0203 function [model1,reversedRxns1] = flipRxn(model1, rxnInd, groupIds1, reversedRxns1) +0204 model1.S(:,rxnInd) = model1.S(:,rxnInd)*-1; +0205 %swap the bounds +0206 ub = model1.ub(rxnInd); +0207 model1.ub(rxnInd) = -model1.lb(rxnInd); +0208 model1.lb(rxnInd) = -ub; +0209 %flip the objective +0210 model1.c(rxnInd) = -model1.c(rxnInd); +0211 +0212 %now take care of the reversedRxns - if this is a group, reverse all of the +0213 %reactions in the group in the reversedRxns index - they will all be reversed at the +0214 %same time since they are the same rxn. +0215 rxnIndices = rxnInd; +0216 if groupIds1(rxnInd) > 0 +0217 rxnIndices = find(groupIds1 == groupIds1(rxnInd)); +0218 end +0219 reversedRxns1(rxnIndices) = ~reversedRxns1(rxnIndices); +0220 end +0221 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 From 0a1b1aafcee0fd880d4b0ba6c596bd52fd5669a8 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 7 Oct 2026 21:07:33 +0200 Subject: [PATCH 3/3] fix: compare rev rather than coerce it, and cover ftINIT too The earlier fix coerced reducedModel.rev to double on entry to mergeLinear, which changed the class of rev in the model it returns -- a type the caller did not ask to have changed. Comparing instead matches what every other use of rev in RAVEN already does (rev == 1, rev ~= 0, find(rev)), works for a double, an integer or a logical alike, and leaves the model untouched. On the released Human-GEM 2.1.0 .mat, whose rev is int64, both forms return the same 8751 reactions with the same S and the same rev values, but the comparison leaves rev as int64 instead of silently making it double. ftINIT had the same pattern and the same failure: any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' a sparse logical combined with rev, which MATLAB refuses for an integer array. It only runs when metabolomics data is supplied, which is why #1143 did not reach it. Every other use of rev compares first and is already class-agnostic. --- INIT/ftINIT.m | 2 +- INIT/mergeLinear.m | 14 +- doc/INIT/ftINIT.html | 2 +- doc/INIT/mergeLinear.html | 356 +++++++++++++++++++------------------- 4 files changed, 185 insertions(+), 189 deletions(-) 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 be642369f..114d7228d 100644 --- a/INIT/mergeLinear.m +++ b/INIT/mergeLinear.m @@ -39,12 +39,6 @@ reducedModel.geneComps=[]; end -%rev is combined with a sparse logical below ('S(...).' > 0 | rev'), which -%MATLAB refuses for an integer class: "Sparse integer array arithmetic -%operations are not supported". A model built outside RAVEN can carry an -%integer rev, so coerce it once here rather than failing several steps in. -reducedModel.rev = double(reducedModel.rev); - nextGroupId = 1; origRxnIds = reducedModel.rxns; groupIds = zeros(numel(reducedModel.rxns),1); @@ -68,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/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 2c5c663a5..5d90456af 100644 --- a/doc/INIT/mergeLinear.html +++ b/doc/INIT/mergeLinear.html @@ -102,186 +102,184 @@

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