diff --git a/INIT/ftINIT.m b/INIT/ftINIT.m index ab5f3e030..9dba76ad1 100644 --- a/INIT/ftINIT.m +++ b/INIT/ftINIT.m @@ -135,7 +135,7 @@ %Get the matching mets metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i})); prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers - (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers + (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' ~= 0); %reversible reactions that are consumers %convert the production rxns from refModel to minModel prepData.groupIds; [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns); diff --git a/INIT/mergeLinear.m b/INIT/mergeLinear.m index 1378353b2..114d7228d 100644 --- a/INIT/mergeLinear.m +++ b/INIT/mergeLinear.m @@ -62,8 +62,12 @@ for i=1:numel(twoNonZero) involvedRxns=find(reducedModel.S(twoNonZero(i),:)); %Check that we can have one positive and one negative - pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns)); - neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns)); + %rev is compared rather than used as-is: the other operand is a + %sparse logical, and MATLAB refuses to combine one with an integer + %array ("Sparse integer array arithmetic operations are not + %supported"). A model built outside RAVEN can carry an integer rev. + pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0); + neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0); %Check so that one or both of the reactions haven't been merged diff --git a/doc/INIT/ftINIT.html b/doc/INIT/ftINIT.html index 6133e894d..8e8a54983 100644 --- a/doc/INIT/ftINIT.html +++ b/doc/INIT/ftINIT.html @@ -270,7 +270,7 @@

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

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

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

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

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

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

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