diff --git a/doc/external/kegg/getPhylDist.html b/doc/external/kegg/getPhylDist.html index 5bf62f6b..8c63971e 100644 --- a/doc/external/kegg/getPhylDist.html +++ b/doc/external/kegg/getPhylDist.html @@ -111,96 +111,98 @@

SOURCE CODE ^'taxonomy'), 'r'); 0049 0050 phylDistStruct.ids={}; -0051 -0052 %Keeps the categories for each organism -0053 orgCat={}; -0054 -0055 currentCat={}; -0056 %Keeps track of the current category -0057 -0058 depth=0; -0059 %Keeps track of the current tree depth -0060 -0061 %Loop through the file -0062 orgCounter=0; -0063 while 1 -0064 %Get the next line -0065 tline = fgetl(fid); -0066 -0067 %Abort at end of file -0068 if ~ischar(tline) -0069 break; -0070 end -0071 -0072 if any(tline) -0073 %Check if it's a new category -0074 if tline(1)=='#' -0075 %Find the first space (=depth +1) -0076 sPos=strfind(tline,' ')-1; -0077 %Should always exist -0078 -0079 sPos=sPos(1); -0080 -0081 %If we have stepped back one step in the tree -0082 if sPos<depth -0083 currentCat=currentCat(1:sPos); -0084 end -0085 depth=sPos; -0086 -0087 currentCat{depth}=tline(sPos+2:end); -0088 else -0089 orgCounter=orgCounter+1; -0090 %It is an organism Get the id between first and second -0091 %white space -0092 sPos=find(isstrprop(tline, 'wspace')); -0093 %Should always exist -0094 -0095 phylDistStruct.ids{orgCounter}=tline(sPos(1)+1:sPos(2)-1); -0096 orgCat{orgCounter}=currentCat; -0097 end -0098 end -0099 end -0100 %Generate a distance matrix (very straight forward here, not neat) -0101 phylDistStruct.distMat=zeros(numel(phylDistStruct.ids)); -0102 phylDistStructOnlyInKingdom.distMat=zeros(numel(phylDistStruct.ids)); -0103 phylDistStructOnlyInKingdom.ids=phylDistStruct.ids; -0104 for i=1:numel(phylDistStruct.ids) -0105 for j=1:numel(phylDistStruct.ids) -0106 if ~strcmp(orgCat{i}(1),orgCat{j}(1)) -0107 phylDistStructOnlyInKingdom.distMat(i,j)=Inf; -0108 end -0109 %Calculate the distance between then -0110 dist=numel(orgCat{i})-numel(orgCat{j}); -0111 if dist>0 -0112 aCat=orgCat{i}(1:end-dist); -0113 else -0114 aCat=orgCat{i}; -0115 end -0116 if dist<0 -0117 bCat=orgCat{j}(1:end+dist); -0118 else -0119 bCat=orgCat{j}; -0120 end -0121 -0122 %Loop through the categories and stop when they are the -0123 %same -0124 for k=numel(aCat):-1:1 -0125 if strcmp(aCat{k},bCat{k}) -0126 break; -0127 end -0128 end -0129 phylDistStruct.distMat(i,j)=dist+numel(aCat)-k; -0130 end -0131 end -0132 %Save the structure -0133 save(distFile,'phylDistStruct','phylDistStructOnlyInKingdom'); -0134 fprintf('COMPLETE\n'); -0135 end -0136 end -0137 if onlyInKingdom==true -0138 phylDistStruct=phylDistStructOnlyInKingdom; -0139 end -0140 end +0051 phylDistStruct.names={}; +0052 +0053 %Keeps the categories for each organism +0054 orgCat={}; +0055 +0056 currentCat={}; +0057 %Keeps track of the current category +0058 +0059 depth=0; +0060 %Keeps track of the current tree depth +0061 +0062 %Loop through the file +0063 orgCounter=0; +0064 while 1 +0065 %Get the next line +0066 tline = fgetl(fid); +0067 +0068 %Abort at end of file +0069 if ~ischar(tline) +0070 break; +0071 end +0072 +0073 if any(tline) +0074 %Check if it's a new category +0075 if tline(1)=='#' +0076 %Find the first space (=depth +1) +0077 sPos=strfind(tline,' ')-1; +0078 %Should always exist +0079 +0080 sPos=sPos(1); +0081 +0082 %If we have stepped back one step in the tree +0083 if sPos<depth +0084 currentCat=currentCat(1:sPos); +0085 end +0086 depth=sPos; +0087 +0088 currentCat{depth}=tline(sPos+2:end); +0089 else +0090 orgCounter=orgCounter+1; +0091 %It is an organism Get the id between first and second +0092 %white space +0093 sPos=find(isstrprop(tline, 'wspace')); +0094 %Should always exist +0095 +0096 phylDistStruct.ids{orgCounter}=tline(sPos(1)+1:sPos(2)-1); +0097 phylDistStruct.names{orgCounter}=tline(sPos(3)+1:end); +0098 orgCat{orgCounter}=currentCat; +0099 end +0100 end +0101 end +0102 %Generate a distance matrix (very straight forward here, not neat) +0103 phylDistStruct.distMat=zeros(numel(phylDistStruct.ids)); +0104 phylDistStructOnlyInKingdom.distMat=zeros(numel(phylDistStruct.ids)); +0105 phylDistStructOnlyInKingdom.ids=phylDistStruct.ids; +0106 for i=1:numel(phylDistStruct.ids) +0107 for j=1:numel(phylDistStruct.ids) +0108 if ~strcmp(orgCat{i}(1),orgCat{j}(1)) +0109 phylDistStructOnlyInKingdom.distMat(i,j)=Inf; +0110 end +0111 %Calculate the distance between then +0112 dist=numel(orgCat{i})-numel(orgCat{j}); +0113 if dist>0 +0114 aCat=orgCat{i}(1:end-dist); +0115 else +0116 aCat=orgCat{i}; +0117 end +0118 if dist<0 +0119 bCat=orgCat{j}(1:end+dist); +0120 else +0121 bCat=orgCat{j}; +0122 end +0123 +0124 %Loop through the categories and stop when they are the +0125 %same +0126 for k=numel(aCat):-1:1 +0127 if strcmp(aCat{k},bCat{k}) +0128 break; +0129 end +0130 end +0131 phylDistStruct.distMat(i,j)=dist+numel(aCat)-k; +0132 end +0133 end +0134 %Save the structure +0135 save(distFile,'phylDistStruct','phylDistStructOnlyInKingdom'); +0136 fprintf('COMPLETE\n'); +0137 end +0138 end +0139 if onlyInKingdom==true +0140 phylDistStruct=phylDistStructOnlyInKingdom; +0141 end +0142 end
Generated by m2html © 2005
\ No newline at end of file diff --git a/doc/external/metacyc/getMetaCycModelForOrganism.html b/doc/external/metacyc/getMetaCycModelForOrganism.html index c664585d..5745903d 100644 --- a/doc/external/metacyc/getMetaCycModelForOrganism.html +++ b/doc/external/metacyc/getMetaCycModelForOrganism.html @@ -24,7 +24,7 @@

PURPOSE ^getMetaCycModelForOrganism

SYNOPSIS ^

-
function model=getMetaCycModelForOrganism(organismID,fastaFile,keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond)
+
function [model, metacycBlastStruct] = getMetaCycModelForOrganism(organismID,fastaFile,keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond,metacycBlastStruct)

DESCRIPTION ^

 getMetaCycModelForOrganism
@@ -50,12 +50,19 @@ 
DESCRIPTION ^CROSS-REFERENCE INFORMATION ^

@@ -70,8 +77,8 @@ 
CROSS-REFERENCE INFORMATION ^
 
 
 <h2><a name=SOURCE CODE ^

-
0001 function model=getMetaCycModelForOrganism(organismID,fastaFile,...
-0002     keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond)
+
0001 function [model, metacycBlastStruct] = getMetaCycModelForOrganism(organismID,fastaFile,...
+0002     keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond,metacycBlastStruct)
 0003 % getMetaCycModelForOrganism
 0004 %   Reconstructs a genome-scale metabolic model based on protein homology to the
 0005 %   MetaCyc pathway database
@@ -95,235 +102,249 @@ 
SOURCE CODE ^%   minPositives        minimum Positives values of BLASTp search (opt, default 45 %)
 0024 %   useDiamond          use DIAMOND alignment tools to perform homology search
 0025 %                       if true, otherwise the BLASTP is used (opt, default true)
-0026 %
-0027 %   Output:
-0028 %   model               a model structure for the query organism
-0029 %
-0030 %   Usage: model=getMetaCycModelForOrganism(organismID,fastaFile,...
-0031 %    keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond)
-0032 
-0033 organismID=char(organismID);
-0034 if nargin<2
-0035     EM='No query protein fasta file is specified';
-0036     dispEM(EM);
-0037 else
-0038     fastaFile=char(fastaFile);
-0039 end
-0040 if nargin<3
-0041     keepTransportRxns=false;
-0042 end
-0043 if nargin<4
-0044     keepUnbalanced=false;
-0045 end
-0046 if nargin<5
-0047     keepUndetermined=false;
-0048 end
-0049 if nargin<6
-0050     minScore=100;
-0051 end
-0052 if nargin<7
-0053     minPositives=45;
-0054 end
-0055 if nargin<8
-0056     useDiamond=true;
-0057 end
-0058 
-0059 if isempty(fastaFile)
-0060     error('*** The query FASTA filename cannot be empty! ***');
-0061 else
-0062     fprintf('\nChecking existence of query FASTA file... ');
-0063     %Check if query fasta exists
-0064     fastaFile=checkFileExistence(fastaFile,2); %Copy file to temp dir
-0065     fprintf('done\n');
-0066 end
-0067 
-0068 %First generate the full MetaCyc model
-0069 metaCycModel=getModelFromMetaCyc([],keepTransportRxns,keepUnbalanced,keepUndetermined);
-0070 fprintf('The full MetaCyc model loaded\n');
-0071 
-0072 %Create the draft model from MetaCyc super model model=metaCycModel;
-0073 model.id=organismID;
-0074 model.name='Generated by homology with MetaCyc database';
-0075 model.rxns=metaCycModel.rxns;
-0076 model.rxnNames=metaCycModel.rxnNames;
-0077 model.eccodes=metaCycModel.eccodes;
-0078 model.subSystems=metaCycModel.subSystems;
-0079 model.rxnMiriams=metaCycModel.rxnMiriams;
-0080 model.rxnReferences=metaCycModel.rxnReferences;
-0081 model.lb=metaCycModel.lb;
-0082 model.ub=metaCycModel.ub;
-0083 model.rev=metaCycModel.rev;
-0084 model.c=metaCycModel.c;
-0085 model.equations=metaCycModel.equations;
-0086 
-0087 %Get the root directory for RAVEN Toolbox.
-0088 ravenPath=findRAVENroot();
-0089 
-0090 %Generate blast strcture by either DIAMOND or BLASTP
-0091 if useDiamond
-0092     blastStruc=getDiamond(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa'));
-0093 else
-0094     blastStruc=getBlast(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa'));
-0095 end
+0026 %   metacycBlastStruct  provided an earlier generated metacycBlastStruct
+0027 %                       from getMetaCycModelForOrganism (opt, default empty).
+0028 %                       Useful when trying different cutoffs for minScore
+0029 %                       and minPositives without having to regenerate the
+0030 %                       blastStructure each time.
+0031 %
+0032 %   Output:
+0033 %   model               a model structure for the query organism
+0034 %   metacycBlastStruct  result from getBlast or getDiamond, before the
+0035 %                       minScore and minPositives cutoffs are applied
+0036 %
+0037 %   Usage: [model, metacycBlastStruct] = getMetaCycModelForOrganism(organismID,fastaFile,...
+0038 %    keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond,metacycBlastStruct)
+0039 
+0040 organismID=char(organismID);
+0041 if nargin<2
+0042     EM='No query protein fasta file is specified';
+0043     dispEM(EM);
+0044 else
+0045     fastaFile=char(fastaFile);
+0046 end
+0047 if nargin<3
+0048     keepTransportRxns=false;
+0049 end
+0050 if nargin<4
+0051     keepUnbalanced=false;
+0052 end
+0053 if nargin<5
+0054     keepUndetermined=false;
+0055 end
+0056 if nargin<6
+0057     minScore=100;
+0058 end
+0059 if nargin<7
+0060     minPositives=45;
+0061 end
+0062 if nargin<8
+0063     useDiamond=true;
+0064 end
+0065 if nargin<9
+0066     metacycBlastStruct=[];
+0067 end
+0068 
+0069 if isempty(fastaFile)
+0070     error('*** The query FASTA filename cannot be empty! ***');
+0071 else
+0072     fprintf('\nChecking existence of query FASTA file... ');
+0073     %Check if query fasta exists
+0074     fastaFile=checkFileExistence(fastaFile,2); %Copy file to temp dir
+0075     fprintf('done\n');
+0076 end
+0077 
+0078 %First generate the full MetaCyc model
+0079 metaCycModel=getModelFromMetaCyc([],keepTransportRxns,keepUnbalanced,keepUndetermined);
+0080 fprintf('The full MetaCyc model loaded\n');
+0081 
+0082 %Create the draft model from MetaCyc super model model=metaCycModel;
+0083 model.id=organismID;
+0084 model.name='Generated by homology with MetaCyc database';
+0085 model.rxns=metaCycModel.rxns;
+0086 model.rxnNames=metaCycModel.rxnNames;
+0087 model.eccodes=metaCycModel.eccodes;
+0088 model.subSystems=metaCycModel.subSystems;
+0089 model.rxnMiriams=metaCycModel.rxnMiriams;
+0090 model.rxnReferences=metaCycModel.rxnReferences;
+0091 model.lb=metaCycModel.lb;
+0092 model.ub=metaCycModel.ub;
+0093 model.rev=metaCycModel.rev;
+0094 model.c=metaCycModel.c;
+0095 model.equations=metaCycModel.equations;
 0096 
-0097 %Only look the query
-0098 blastStructure=blastStruc(2);
+0097 %Get the root directory for RAVEN Toolbox.
+0098 ravenPath=findRAVENroot();
 0099 
-0100 %Remove all blast hits that are below the cutoffs
-0101 indexes=blastStructure.bitscore>=minScore & blastStructure.ppos>=minPositives;
-0102 blastStructure.toGenes(~indexes)=[];
-0103 blastStructure.fromGenes(~indexes)=[];
-0104 blastStructure.evalue(~indexes)=[];
-0105 blastStructure.aligLen(~indexes)=[];
-0106 blastStructure.identity(~indexes)=[];
-0107 blastStructure.bitscore(~indexes)=[];
-0108 blastStructure.ppos(~indexes)=[];
-0109 fprintf('Completed searching against MetaCyc protein sequences.\n');
-0110 
-0111 % Get the qualified genes of query organism from blast structure
-0112 model.genes=cell(100000,1);
-0113 model.proteins=cell(100000,1);
-0114 model.bitscore=zeros(100000,1);
-0115 model.ppos=zeros(100000,1);
-0116 num=1;
-0117 
-0118 %Go through the strucutre and find out the hit with the best bit score
-0119 for i=1:numel(blastStructure.fromGenes)
-0120     if num==1
-0121         model.genes(num)=blastStructure.fromGenes(i);
-0122         model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
-0123         model.bitscore(num)=blastStructure.bitscore(i);
-0124         model.ppos(num)=blastStructure.ppos(i);
-0125         num=num+1;
-0126         lastGene=blastStructure.fromGenes(i);
-0127     else
-0128         if ~isequal(lastGene,blastStructure.fromGenes(i))
-0129             model.genes(num)=blastStructure.fromGenes(i);
-0130             model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
-0131             model.bitscore(num)=blastStructure.bitscore(i);
-0132             model.ppos(num)=blastStructure.ppos(i);
-0133             num=num+1;
-0134             lastGene=blastStructure.fromGenes(i);
-0135         else
-0136             if model.bitscore(num)<blastStructure.bitscore(i)
-0137                 model.bitscore(num)=blastStructure.bitscore(i);
-0138                 model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
-0139                 model.ppos(num)=blastStructure.ppos(i);
-0140             end
-0141         end
-0142     end
-0143 end
-0144 model.genes=model.genes(1:num-1);
-0145 model.proteins=model.proteins(1:num-1);
-0146 model.bitscore=model.bitscore(1:num-1);
-0147 model.ppos=model.ppos(1:num-1);
-0148 
-0149 % Get the indexes of matched genes in the metaCycModel
-0150 % because some enzymes in the FASTA file cannot be found in the dump file
-0151 [hits, indexes]=ismember(model.proteins,metaCycModel.genes);
-0152 found = find(hits);
-0153 model.genes=model.genes(found);
-0154 
-0155 % Restructure the rxnGeneMat matrix
-0156 model.rxnGeneMat=sparse(metaCycModel.rxnGeneMat(:,indexes(found)));
-0157 
-0158 %Remove all reactions without genes
-0159 hasGenes=any(model.rxnGeneMat,2);
-0160 model=removeReactions(model,~hasGenes,true);
-0161 
-0162 %Generate grRules, only consider the or relationship here Matched enzymes
-0163 %are stored in field rxnScores,
-0164 rxnNum=numel(model.rxns);
-0165 model.rxnConfidenceScores=NaN(rxnNum,1);
-0166 model.rxnConfidenceScores(:)=2;
-0167 model.grRules=cell(rxnNum,1);
-0168 %model.rxnScores=cell(rxnNum,1);
-0169 for j=1:rxnNum
-0170     model.grRules{j}='';
-0171     I=find(model.rxnGeneMat(j,:));
-0172     for k=1:numel(I)
-0173         if isempty(model.grRules{j})
-0174             model.grRules(j)=model.genes(I(k));
-0175             %model.rxnScores(j)=model.proteins(I(k));
-0176         else
-0177             model.grRules(j)=strcat(model.grRules(j),{' or '},model.genes(I(k)));
-0178             %model.rxnScores(j)=strcat(model.rxnScores(j),{' or
-0179             %'},model.proteins(I(k)));
-0180         end
-0181     end
-0182 end
-0183 %update genes field
-0184 model.genes=model.genes(any(model.rxnGeneMat));
-0185 
-0186 %Construct the S matrix and list of metabolites
-0187 [S, mets, badRxns]=constructS(model.equations);
-0188 model.S=S;
-0189 model.mets=mets;
-0190 
-0191 %model=removeReactions(model,badRxns,true,true);
-0192 
-0193 %Then match up metabolites
-0194 metaCycMets=getMetsFromMetaCyc([]);
-0195 
-0196 %Add information about all metabolites to the model
-0197 [a, b]=ismember(model.mets,metaCycMets.mets);
-0198 a=find(a);
-0199 b=b(a);
-0200 
-0201 if ~isfield(model,'metNames')
-0202     model.metNames=cell(numel(model.mets),1);
-0203     model.metNames(:)={''};
-0204 end
-0205 model.metNames(a)=metaCycMets.metNames(b);
+0100 %Generate blast strcture by either DIAMOND or BLASTP
+0101 if isempty(metacycBlastStruct)
+0102     if useDiamond
+0103         blastStruc=getDiamond(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa'));
+0104     else
+0105         blastStruc=getBlast(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa'));
+0106     end
+0107     %Only look the query
+0108     blastStructure=blastStruc(2);
+0109     metacycBlastStruct=blastStructure;
+0110 else
+0111     blastStructure=metacycBlastStruct;
+0112 end
+0113 
+0114 %Remove all blast hits that are below the cutoffs
+0115 indexes=blastStructure.bitscore>=minScore & blastStructure.ppos>=minPositives;
+0116 blastStructure.toGenes(~indexes)=[];
+0117 blastStructure.fromGenes(~indexes)=[];
+0118 blastStructure.evalue(~indexes)=[];
+0119 blastStructure.aligLen(~indexes)=[];
+0120 blastStructure.identity(~indexes)=[];
+0121 blastStructure.bitscore(~indexes)=[];
+0122 blastStructure.ppos(~indexes)=[];
+0123 fprintf('Completed searching against MetaCyc protein sequences.\n');
+0124 
+0125 % Get the qualified genes of query organism from blast structure
+0126 model.genes=cell(100000,1);
+0127 model.proteins=cell(100000,1);
+0128 model.bitscore=zeros(100000,1);
+0129 model.ppos=zeros(100000,1);
+0130 num=1;
+0131 
+0132 %Go through the strucutre and find out the hit with the best bit score
+0133 for i=1:numel(blastStructure.fromGenes)
+0134     if num==1
+0135         model.genes(num)=blastStructure.fromGenes(i);
+0136         model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
+0137         model.bitscore(num)=blastStructure.bitscore(i);
+0138         model.ppos(num)=blastStructure.ppos(i);
+0139         num=num+1;
+0140         lastGene=blastStructure.fromGenes(i);
+0141     else
+0142         if ~isequal(lastGene,blastStructure.fromGenes(i))
+0143             model.genes(num)=blastStructure.fromGenes(i);
+0144             model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
+0145             model.bitscore(num)=blastStructure.bitscore(i);
+0146             model.ppos(num)=blastStructure.ppos(i);
+0147             num=num+1;
+0148             lastGene=blastStructure.fromGenes(i);
+0149         else
+0150             if model.bitscore(num)<blastStructure.bitscore(i)
+0151                 model.bitscore(num)=blastStructure.bitscore(i);
+0152                 model.proteins(num)=strrep(blastStructure.toGenes(i), 'gnl|META|', '');
+0153                 model.ppos(num)=blastStructure.ppos(i);
+0154             end
+0155         end
+0156     end
+0157 end
+0158 model.genes=model.genes(1:num-1);
+0159 model.proteins=model.proteins(1:num-1);
+0160 model.bitscore=model.bitscore(1:num-1);
+0161 model.ppos=model.ppos(1:num-1);
+0162 
+0163 % Get the indexes of matched genes in the metaCycModel
+0164 % because some enzymes in the FASTA file cannot be found in the dump file
+0165 [hits, indexes]=ismember(model.proteins,metaCycModel.genes);
+0166 found = find(hits);
+0167 model.genes=model.genes(found);
+0168 
+0169 % Restructure the rxnGeneMat matrix
+0170 model.rxnGeneMat=sparse(metaCycModel.rxnGeneMat(:,indexes(found)));
+0171 
+0172 %Remove all reactions without genes
+0173 hasGenes=any(model.rxnGeneMat,2);
+0174 model=removeReactions(model,~hasGenes,true);
+0175 
+0176 %Generate grRules, only consider the or relationship here Matched enzymes
+0177 %are stored in field rxnScores,
+0178 rxnNum=numel(model.rxns);
+0179 model.rxnConfidenceScores=NaN(rxnNum,1);
+0180 model.rxnConfidenceScores(:)=2;
+0181 model.grRules=cell(rxnNum,1);
+0182 %model.rxnScores=cell(rxnNum,1);
+0183 for j=1:rxnNum
+0184     model.grRules{j}='';
+0185     I=find(model.rxnGeneMat(j,:));
+0186     for k=1:numel(I)
+0187         if isempty(model.grRules{j})
+0188             model.grRules(j)=model.genes(I(k));
+0189             %model.rxnScores(j)=model.proteins(I(k));
+0190         else
+0191             model.grRules(j)=strcat(model.grRules(j),{' or '},model.genes(I(k)));
+0192             %model.rxnScores(j)=strcat(model.rxnScores(j),{' or
+0193             %'},model.proteins(I(k)));
+0194         end
+0195     end
+0196 end
+0197 %update genes field
+0198 model.genes=model.genes(any(model.rxnGeneMat));
+0199 
+0200 %Construct the S matrix and list of metabolites
+0201 [S, mets, badRxns]=constructS(model.equations);
+0202 model.S=S;
+0203 model.mets=mets;
+0204 
+0205 %model=removeReactions(model,badRxns,true,true);
 0206 
-0207 if ~isfield(model,'metFormulas')
-0208     model.metFormulas=cell(numel(model.mets),1);
-0209     model.metFormulas(:)={''};
-0210 end
-0211 model.metFormulas(a)=metaCycMets.metFormulas(b);
-0212 
-0213 if ~isfield(model,'metCharges')
-0214     model.metCharges=zeros(numel(model.mets),1);
-0215 end
-0216 model.metCharges(a)=metaCycMets.metCharges(b);
-0217 
-0218 if ~isfield(model,'b')
-0219     model.b=zeros(numel(model.mets),1);
-0220 end
-0221 %model.b(a)=metaCycMets.b(b);
-0222 
-0223 if ~isfield(model,'inchis')
-0224     model.inchis=cell(numel(model.mets),1);
-0225     model.inchis(:)={''};
-0226 end
-0227 model.inchis(a)=metaCycMets.inchis(b);
-0228 
-0229 if ~isfield(model,'metMiriams')
-0230     model.metMiriams=cell(numel(model.mets),1);
-0231 end
-0232 model.metMiriams(a)=metaCycMets.metMiriams(b);
-0233 
-0234 %Put all metabolites in one compartment called 's' (for system). This is
-0235 %done just to be more compatible with the rest of the code
-0236 model.comps={'s'};
-0237 model.compNames={'System'};
-0238 model.metComps=ones(numel(model.mets),1);
-0239 
-0240 %It could also be that the metabolite names are empty for some reason In
-0241 %that case, use the ID instead
-0242 I=cellfun(@isempty,model.metNames);
-0243 model.metNames(I)=model.mets(I);
-0244 
-0245 %Remove additional fields
-0246 model=rmfield(model,{'proteins','bitscore','ppos'});
+0207 %Then match up metabolites
+0208 metaCycMets=getMetsFromMetaCyc([]);
+0209 
+0210 %Add information about all metabolites to the model
+0211 [a, b]=ismember(model.mets,metaCycMets.mets);
+0212 a=find(a);
+0213 b=b(a);
+0214 
+0215 if ~isfield(model,'metNames')
+0216     model.metNames=cell(numel(model.mets),1);
+0217     model.metNames(:)={''};
+0218 end
+0219 model.metNames(a)=metaCycMets.metNames(b);
+0220 
+0221 if ~isfield(model,'metFormulas')
+0222     model.metFormulas=cell(numel(model.mets),1);
+0223     model.metFormulas(:)={''};
+0224 end
+0225 model.metFormulas(a)=metaCycMets.metFormulas(b);
+0226 
+0227 if ~isfield(model,'metCharges')
+0228     model.metCharges=zeros(numel(model.mets),1);
+0229 end
+0230 model.metCharges(a)=metaCycMets.metCharges(b);
+0231 
+0232 if ~isfield(model,'b')
+0233     model.b=zeros(numel(model.mets),1);
+0234 end
+0235 %model.b(a)=metaCycMets.b(b);
+0236 
+0237 if ~isfield(model,'inchis')
+0238     model.inchis=cell(numel(model.mets),1);
+0239     model.inchis(:)={''};
+0240 end
+0241 model.inchis(a)=metaCycMets.inchis(b);
+0242 
+0243 if ~isfield(model,'metMiriams')
+0244     model.metMiriams=cell(numel(model.mets),1);
+0245 end
+0246 model.metMiriams(a)=metaCycMets.metMiriams(b);
 0247 
-0248 %In the end fix grRules and rxnGeneMat
-0249 [grRules,rxnGeneMat] = standardizeGrRules(model,false); %Get detailed output
-0250 model.grRules = grRules;
-0251 model.rxnGeneMat = rxnGeneMat;
-0252 %Remove the temp fasta file
-0253 delete(fastaFile)
-0254 end

+0248 %Put all metabolites in one compartment called 's' (for system). This is
+0249 %done just to be more compatible with the rest of the code
+0250 model.comps={'s'};
+0251 model.compNames={'System'};
+0252 model.metComps=ones(numel(model.mets),1);
+0253 
+0254 %It could also be that the metabolite names are empty for some reason In
+0255 %that case, use the ID instead
+0256 I=cellfun(@isempty,model.metNames);
+0257 model.metNames(I)=model.mets(I);
+0258 
+0259 %Remove additional fields
+0260 model=rmfield(model,{'proteins','bitscore','ppos'});
+0261 
+0262 %In the end fix grRules and rxnGeneMat
+0263 [grRules,rxnGeneMat] = standardizeGrRules(model,false); %Get detailed output
+0264 model.grRules = grRules;
+0265 model.rxnGeneMat = rxnGeneMat;
+0266 %Remove the temp fasta file
+0267 delete(fastaFile)
+0268 end

 
Generated by m2html © 2005

 
 
\ No newline at end of file
diff --git a/doc/external/metacyc/getMetsFromMetaCyc.html b/doc/external/metacyc/getMetsFromMetaCyc.html
index 73ddd0b5..97d2a9e9 100644
--- a/doc/external/metacyc/getMetsFromMetaCyc.html
+++ b/doc/external/metacyc/getMetsFromMetaCyc.html
@@ -66,7 +66,7 @@ 
CROSS-REFERENCE INFORMATION ^
 </ul>
 This function is called by:
 <ul style=
-
  • getMetaCycModelForOrganism	getMetaCycModelForOrganism
  • getModelFromMetaCyc	getModelFromMetaCyc
    • 
+
  • getMetaCycModelForOrganism	getMetaCycModelForOrganism
  • getModelFromMetaCyc	getModelFromMetaCyc
    • 
 
 
 
diff --git a/doc/external/metacyc/getModelFromMetaCyc.html b/doc/external/metacyc/getModelFromMetaCyc.html
index 93f99558..81a0a226 100644
--- a/doc/external/metacyc/getModelFromMetaCyc.html
+++ b/doc/external/metacyc/getModelFromMetaCyc.html
@@ -68,7 +68,7 @@ 
    CROSS-REFERENCE INFORMATION ^
 <li><a href=getEnzymesFromMetaCyc	getEnzymesFromMetaCyc
  • getMetsFromMetaCyc	getMetsFromMetaCyc
  • getRxnsFromMetaCyc	getRxnsFromMetaCyc
    • 
 This function is called by:
 
+
  • getMetaCycModelForOrganism	getMetaCycModelForOrganism
    • 
 
 
 
diff --git a/doc/installation/checkFunctionUniqueness.html b/doc/installation/checkFunctionUniqueness.html
index ad22f841..a2075ffe 100644
--- a/doc/installation/checkFunctionUniqueness.html
+++ b/doc/installation/checkFunctionUniqueness.html
@@ -79,78 +79,83 @@ 
    SOURCE CODE ^else
-0026     altDirs = strsplit(altDirs,';')';
-0027     ravenFunctions={};
-0028     for i=1:numel(altDirs)
-0029         temp_res=dir(fullfile(char(altDirs(i)),'*.m'));
-0030         ravenFunctions = [ravenFunctions, temp_res.name];
-0031     end
-0032     ravenFunctions=ravenFunctions';
-0033     ravenDir=altDirs;
-0034 end
-0035 
-0036 %startup.m is not a normal function, any startup.m in the path should run
-0037 %during startup, so duplicate use of this name is fine
-0038 ravenFunctions=ravenFunctions(~ismember(ravenFunctions,'startup.m'));
-0039 
-0040 %Getting all the paths added to Matlab
-0041 if ispc
-0042     matlabPaths=regexp(path, ';', 'split')';
-0043 elseif isunix
-0044     matlabPaths=regexp(path, ':', 'split')';
-0045 end
-0046 
-0047 overlapPath={};
-0048 overlapFunctions={};
-0049 multiRaven=false;
-0050 multiFunction=false;
+0026     if ispc
+0027         splitChar = ';';
+0028     else
+0029         splitChar = ':';
+0030     end
+0031     altDirs = strsplit(altDirs,splitChar)';
+0032     ravenFunctions={};
+0033     for i=1:numel(altDirs)
+0034         temp_res=dir(fullfile(char(altDirs(i)),'*.m'));
+0035         ravenFunctions = [ravenFunctions, temp_res.name];
+0036     end
+0037     ravenFunctions=ravenFunctions';
+0038     ravenDir=altDirs;
+0039 end
+0040 
+0041 %startup.m is not a normal function, any startup.m in the path should run
+0042 %during startup, so duplicate use of this name is fine
+0043 ravenFunctions=ravenFunctions(~ismember(ravenFunctions,'startup.m'));
+0044 
+0045 %Getting all the paths added to Matlab
+0046 if ispc
+0047     matlabPaths=regexp(path, ';', 'split')';
+0048 elseif isunix
+0049     matlabPaths=regexp(path, ':', 'split')';
+0050 end
 0051 
-0052 for i=1:numel(matlabPaths)
-0053     if ~startsWith(matlabPaths{i},ravenDir)
-0054         temp_res=dir([matlabPaths{i} '/*.m']);
-0055         if ~isempty(temp_res)
-0056             pathFunctions={temp_res.name}';
-0057         else
-0058             pathFunctions='';
-0059         end
-0060         if ~isempty(pathFunctions) && ~any(ismember('Contents.m',pathFunctions))
-0061             if any(ismember(ravenFunctions,pathFunctions))
-0062                 if sum(ismember(ravenFunctions,pathFunctions))>(numel(ravenFunctions)/4)
-0063                     multiRaven=true;
-0064                 else
-0065                     multiFunction=true;
-0066                     overlapPath = [overlapPath, matlabPaths(i)];
-0067                     overlapFunctions = [overlapFunctions, {ravenFunctions(ismember(ravenFunctions,pathFunctions))}];
-0068                 end
-0069             end
-0070         end
-0071     end
-0072 end
-0073 
-0074 if multiRaven==true || multiFunction == true
-0075         if nargout > 0
-0076             status = false;
-0077         else
-0078             fprintf('Fail\n')
-0079         end
-0080     if multiRaven==true && nargin < 1
-0081         error(['Multiple RAVEN versions detected in MATLAB path. Remove all ',...
-0082                'RAVEN directories from the MATLAB path with removeRavenFromPath(), ',...
-0083                'or manually remove them. Afterwards, re-run checkInstallation']);
-0084     elseif multiFunction == true
-0085         for i=1:numel(overlapPath)
-0086             fprintf(['   Duplicate functions in ',regexprep(overlapPath{i},'(\\)','\\$1'),'\n']);
-0087             fprintf(['     ' strjoin(overlapFunctions{i},'\n     ') '\n']);
-0088         end
-0089     fprintf('   Resolve conflicting functions to avoid unexpected behaviour\n');
-0090     end
-0091 else
-0092     if nargout > 0
-0093         status = true;
-0094     else
-0095         fprintf('Pass\n')
-0096     end
-0097 end
    +0052 overlapPath={}; +0053 overlapFunctions={}; +0054 multiRaven=false; +0055 multiFunction=false; +0056 +0057 for i=1:numel(matlabPaths) +0058 if ~startsWith(matlabPaths{i},ravenDir) +0059 temp_res=dir([matlabPaths{i} '/*.m']); +0060 if ~isempty(temp_res) +0061 pathFunctions={temp_res.name}'; +0062 else +0063 pathFunctions=''; +0064 end +0065 if ~isempty(pathFunctions) && ~any(ismember('Contents.m',pathFunctions)) +0066 if any(ismember(ravenFunctions,pathFunctions)) +0067 if sum(ismember(ravenFunctions,pathFunctions))>(numel(ravenFunctions)/4) +0068 multiRaven=true; +0069 else +0070 multiFunction=true; +0071 overlapPath = [overlapPath, matlabPaths(i)]; +0072 overlapFunctions = [overlapFunctions, {ravenFunctions(ismember(ravenFunctions,pathFunctions))}]; +0073 end +0074 end +0075 end +0076 end +0077 end +0078 +0079 if multiRaven==true || multiFunction == true +0080 if nargout > 0 +0081 status = false; +0082 else +0083 fprintf('Fail\n') +0084 end +0085 if multiRaven==true && nargin < 1 +0086 error(['Multiple RAVEN versions detected in MATLAB path. Remove all ',... +0087 'RAVEN directories from the MATLAB path with removeRavenFromPath(), ',... +0088 'or manually remove them. Afterwards, re-run checkInstallation']); +0089 elseif multiFunction == true +0090 for i=1:numel(overlapPath) +0091 fprintf([' Duplicate functions in ',regexprep(overlapPath{i},'(\\)','\\$1'),'\n']); +0092 fprintf([' ' strjoin(overlapFunctions{i},'\n ') '\n']); +0093 end +0094 fprintf(' Resolve conflicting functions to avoid unexpected behaviour\n'); +0095 end +0096 else +0097 if nargout > 0 +0098 status = true; +0099 else +0100 fprintf('Pass\n') +0101 end +0102 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/installation/checkInstallation.html b/doc/installation/checkInstallation.html index 14b0d50e..a9336bb5 100644 --- a/doc/installation/checkInstallation.html +++ b/doc/installation/checkInstallation.html @@ -57,7 +57,7 @@

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

    +
  • function res = interpretResults(results)
  • function str = myStr(InputStr,len)
  • function status = makeBinaryExecutable()
  • function printOrange(stringToPrint)

    • SOURCE CODE ^

    0001 function currVer = checkInstallation(developMode)
@@ -112,13 +112,13 @@ 
    SOURCE CODE ^for i=1:3
 0051                 if currVerNum(i)<newVerNum(i)
 0052                     fprintf([myStr('   > Latest RAVEN release available',40) '%f'])
-0053                     printOrange([newVer,'\n'])
+0053                     printOrange([newVer,'\n'])
 0054                     hasGit=isfolder(fullfile(ravenDir,'.git'));
 0055                     if hasGit
-0056                         printOrange('     Run git pull in your favourite git client\n')
-0057                         printOrange('     to get the latest RAVEN release\n');
+0056                         printOrange('     Run git pull in your favourite git client\n')
+0057                         printOrange('     to get the latest RAVEN release\n');
 0058                     else
-0059                         printOrange([myStr('     Instructions on how to upgrade',40) '%f'])
+0059                         printOrange([myStr('     Instructions on how to upgrade',40) '%f'])
 0060                         fprintf('<a href="https://github.com/SysBioChalmers/RAVEN/wiki/Installation#upgrade-to-latest-raven-release">here</a>\n');
 0061                     end
 0062                     break
@@ -128,8 +128,8 @@ 
    SOURCE CODE ^end
 0067         catch
 0068             fprintf([myStr('   > Checking for latest RAVEN release',40) '%f'])
-0069             printOrange('Fail\n');
-0070             printOrange('     Cannot reach GitHub for release info\n');
+0069             printOrange('Fail\n');
+0070             printOrange('     Cannot reach GitHub for release info\n');
 0071         end
 0072     end
 0073 else
@@ -156,12 +156,12 @@ 
    SOURCE CODE ^'Pass\n')   
 0096     catch
-0097         printOrange('Fail\n')
+0097         printOrange('Fail\n')
 0098         fprintf(['   You might have to rerun checkInstallation again\n'...
 0099                  '   next time you start up MATLAB\n'])        
 0100     end
 0101 catch
-0102     printOrange('Fail\n')
+0102     printOrange('Fail\n')
 0103 end
 0104 
 0105 if isunix
@@ -170,7 +170,7 @@ 
    SOURCE CODE ^if status == 0
 0109         fprintf('Pass\n')
 0110     else
-0111         printOrange('Fail\n')
+0111         printOrange('Fail\n')
 0112     end
 0113 end
 0114 
@@ -182,7 +182,7 @@ 
    SOURCE CODE ^'Pass\n')
 0122 catch
-0123     printOrange('Fail\n')
+0123     printOrange('Fail\n')
 0124 end
 0125 fprintf([myStr(' > Checking libSBML version',40) '%f'])
 0126 try
@@ -191,11 +191,11 @@ 
    SOURCE CODE ^% Only works in libSBML 5.17.0+
 0130         fprintf([libSBMLver.libSBML_version_string '\n']);
 0131     catch
-0132         printOrange('Fail\n')
+0132         printOrange('Fail\n')
 0133         fprintf('   An older libSBML version was found, update to version 5.17.0 or higher for a significant improvement of model import\n');
 0134     end
 0135 catch
-0136     printOrange('Fail\n')
+0136     printOrange('Fail\n')
 0137     fprintf('   Download libSBML from http://sbml.org/Software/libSBML/Downloading_libSBML and add to MATLAB path\n');
 0138 end
 0139 fprintf(' > Checking model import and export\n')
@@ -205,7 +205,7 @@ 
    SOURCE CODE ^if res(1).Passed == 1
 0144     fprintf('Pass\n')
 0145 else
-0146     printOrange('Fail\n')
+0146     printOrange('Fail\n')
 0147     addList = matlab.addons.installedAddons;
 0148     if any(strcmpi(addList.Name,'Text Analytics Toolbox'))
 0149         fprintf(['   Excel import/export is incompatible with MATLAB Text Analytics Toolbox.\n' ...
@@ -217,21 +217,21 @@ 
    SOURCE CODE ^if res(3).Passed == 1
 0156     fprintf('Pass\n')
 0157 else
-0158     printOrange('Fail\n')
+0158     printOrange('Fail\n')
 0159 end
 0160 
 0161 fprintf([myStr('   > Import SBML format',40) '%f'])
 0162 if res(2).Passed == 1
 0163     fprintf('Pass\n')
 0164 else
-0165     printOrange('Fail\n')
+0165     printOrange('Fail\n')
 0166 end
 0167 
 0168 fprintf([myStr('   > Export SBML format',40) '%f'])
 0169 if res(4).Passed == 1
 0170     fprintf('Pass\n')
 0171 else
-0172     printOrange('Fail\n')
+0172     printOrange('Fail\n')
 0173 end
 0174 
 0175 if res(1).Passed~=1 && res(3).Passed~=1 && exist('vaderSentimentScores.m','file')==2
@@ -259,28 +259,28 @@ 
    SOURCE CODE ^if res(1).Passed == 1
 0198     fprintf('Pass\n')
 0199 else
-0200     printOrange('Fail\n')
+0200     printOrange('Fail\n')
 0201 end
 0202 
 0203 fprintf([myStr('   > gurobi',40) '%f'])
 0204 if res(2).Passed == 1
 0205     fprintf('Pass\n')
 0206 else
-0207     printOrange('Fail\n')
+0207     printOrange('Fail\n')
 0208 end
 0209 
 0210 fprintf([myStr('   > soplex',40) '%f'])
 0211 if res(3).Passed == 1
 0212     fprintf('Pass\n')
 0213 else
-0214     printOrange('Fail\n')
+0214     printOrange('Fail\n')
 0215 end
 0216 
 0217 fprintf([myStr('   > cobra',40) '%f'])
 0218 if res(4).Passed == 1
 0219     fprintf('Pass\n')
 0220 else
-0221     printOrange('Fail\n')
+0221     printOrange('Fail\n')
 0222 end
 0223 fprintf([myStr(' > Set RAVEN solver',40) '%f'])
 0224 try
@@ -358,7 +358,7 @@ 
    SOURCE CODE ^'Pass\n');
 0297     res=true;
 0298 else
-0299     printOrange('Fail\n')
+0299     printOrange('Fail\n')
 0300     fprintf('   Download/compile the binary and rerun checkInstallation\n');
 0301     res=false;
 0302 end
@@ -400,7 +400,19 @@ 
    SOURCE CODE ^'Failed to make %s executable: %s ',binList{i},strip(cmdout))
 0339     end
 0340 end
-0341 end
    +0341 end +0342 +0343 function printOrange(stringToPrint) +0344 % printOrange +0345 % Duplicate of RAVEN/core/printOrange is also kept here, as this function +0346 % should be able to run before adding RAVEN to the MATLAB path. +0347 try useDesktop = usejava('desktop'); catch, useDesktop = false; end +0348 if useDesktop +0349 fprintf(['[\b' stringToPrint,']\b']) +0350 else +0351 fprintf(stringToPrint) +0352 end +0353 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/importModel.html b/doc/io/importModel.html index 7e072861..63e03584 100644 --- a/doc/io/importModel.html +++ b/doc/io/importModel.html @@ -644,699 +644,705 @@

    SOURCE CODE ^end -0528 reactionLB(counter)=str2num(lb); -0529 reactionUB(counter)=str2num(ub); -0530 %The order of these parameters should not be hard coded -0531 elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter') -0532 reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value; -0533 reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value; -0534 reactionObjective(counter)=modelSBML.reaction(i).kineticLaw.parameter(3).value; -0535 else -0536 if reactionReversibility(counter)==true -0537 reactionLB(counter)=-inf; -0538 else -0539 reactionLB(counter)=0; -0540 end -0541 reactionUB(counter)=inf; -0542 reactionObjective(counter)=0; -0543 end -0544 -0545 %Find the associated gene if available -0546 %If FBC, get gene association data from corresponding fields -0547 if isfield(modelSBML.reaction(i),'fbc_geneProductAssociation') -0548 if ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation) && ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association) -0549 grRules{counter}=modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association; -0550 end -0551 elseif isfield(modelSBML.reaction(i),'notes') -0552 %This section was previously executed only if isSBML2COBRA is true. Now -0553 %it will be executed, if 'GENE_ASSOCIATION' is found in -0554 %modelSBML.reaction(i).notes -0555 if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION') -0556 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION'); -0557 elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION') -0558 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION'); -0559 else -0560 geneAssociation=''; -0561 end -0562 if ~isempty(geneAssociation) -0563 %This adds the grRules. The gene list and rxnGeneMat are created -0564 %later -0565 grRules{counter}=geneAssociation; -0566 end -0567 end -0568 if isempty(grRules{counter}) && ~isempty(modelSBML.reaction(i).modifier) -0569 rules=''; -0570 for j=1:numel(modelSBML.reaction(i).modifier) -0571 modifier=modelSBML.reaction(i).modifier(j).species; -0572 if ~isempty(modifier) -0573 if strcmpi(modifier(1:2),'E_') -0574 index=find(strcmp(modifier,geneIDs)); -0575 %This should be unique and in the geneIDs list, -0576 %otherwise something is wrong -0577 if numel(index)~=1 -0578 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; -0579 dispEM(EM); -0580 end -0581 if ~isempty(rules) -0582 rules=[rules ' or (' geneNames{index} ')']; -0583 else -0584 rules=['(' geneNames{index} ')']; -0585 end -0586 elseif strcmp(modifier(1:2),'s_') -0587 index=find(strcmp(modifier,metaboliteIDs)); -0588 %This should be unique and in the geneIDs list, -0589 %otherwise something is wrong -0590 if numel(index)~=1 -0591 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; -0592 dispEM(EM); -0593 end -0594 if ~isempty(rules) -0595 rules=[rules ' or (' metaboliteIDs{index} ')']; -0596 else -0597 rules=['(' metaboliteIDs{index} ')']; -0598 end -0599 else -0600 %It seems to be a complex. Add the corresponding -0601 %genes from the name of the complex (not the -0602 %reaction that creates it) -0603 index=find(strcmp(modifier,complexIDs)); -0604 if numel(index)==1 -0605 if ~isempty(rules) -0606 rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')']; -0607 else -0608 rules=['(' strrep(complexNames{index},':',' and ') ')']; -0609 end -0610 else -0611 %Could not find a complex -0612 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; -0613 dispEM(EM); -0614 end -0615 end -0616 end -0617 end -0618 grRules{counter}=rules; -0619 grRulesFromModifier{counter}=rules;%Backup copy for grRules, useful to parse Yeast 7.6 -0620 end -0621 -0622 %Add reaction compartment -0623 if isfield(modelSBML.reaction(i),'compartment') -0624 if ~isempty(modelSBML.reaction(i).compartment) -0625 rxnComp=modelSBML.reaction(i).compartment; -0626 else -0627 rxnComp=''; -0628 end -0629 elseif isfield(modelSBML.reaction(i),'notes') -0630 rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT'); -0631 end -0632 if ~isempty(rxnComp) -0633 %Find it in the compartment list -0634 [~, J]=ismember(rxnComp,compartmentIDs); -0635 rxnComps(counter)=J; -0636 end -0637 -0638 %Get other Miriam fields. This may include for example database indexes -0639 %to organism-specific databases. EC-codes are supported by the COBRA -0640 %Toolbox format and are therefore loaded separately -0641 if isSBML2COBRA==false -0642 miriamStruct=parseMiriam(modelSBML.reaction(i).annotation); -0643 rxnMiriams{counter}=miriamStruct; -0644 if isfield(modelSBML.reaction(i),'notes') -0645 subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM')); -0646 subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = []; -0647 if strfind(modelSBML.reaction(i).notes,'Confidence Level') -0648 rxnconfidencescores(counter)=str2num(parseNote(modelSBML.reaction(i).notes,'Confidence Level')); -0649 end -0650 rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS'); -0651 rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES'); -0652 end -0653 end -0654 -0655 %Get SBO terms -0656 if isfield(modelSBML.reaction(i),'sboTerm') -0657 rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm); -0658 end -0659 -0660 %Get ec-codes -0661 eccode=''; -0662 if ~isempty(modelSBML.reaction(i).annotation) -0663 if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code') -0664 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code'); -0665 elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code') -0666 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code'); -0667 elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code') -0668 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code'); -0669 end -0670 elseif isfield(modelSBML.reaction(i),'notes') -0671 if strfind(modelSBML.reaction(i).notes,'EC Number') -0672 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')]; -0673 elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS') -0674 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')]; +0528 if isempty(lb) +0529 lb='-Inf'; +0530 end +0531 if isempty(ub) +0532 ub='Inf'; +0533 end +0534 reactionLB(counter)=str2num(lb); +0535 reactionUB(counter)=str2num(ub); +0536 %The order of these parameters should not be hard coded +0537 elseif isfield(modelSBML.reaction(i).kineticLaw,'parameter') +0538 reactionLB(counter)=modelSBML.reaction(i).kineticLaw.parameter(1).value; +0539 reactionUB(counter)=modelSBML.reaction(i).kineticLaw.parameter(2).value; +0540 reactionObjective(counter)=modelSBML.reaction(i).kineticLaw.parameter(3).value; +0541 else +0542 if reactionReversibility(counter)==true +0543 reactionLB(counter)=-inf; +0544 else +0545 reactionLB(counter)=0; +0546 end +0547 reactionUB(counter)=inf; +0548 reactionObjective(counter)=0; +0549 end +0550 +0551 %Find the associated gene if available +0552 %If FBC, get gene association data from corresponding fields +0553 if isfield(modelSBML.reaction(i),'fbc_geneProductAssociation') +0554 if ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation) && ~isempty(modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association) +0555 grRules{counter}=modelSBML.reaction(i).fbc_geneProductAssociation.fbc_association.fbc_association; +0556 end +0557 elseif isfield(modelSBML.reaction(i),'notes') +0558 %This section was previously executed only if isSBML2COBRA is true. Now +0559 %it will be executed, if 'GENE_ASSOCIATION' is found in +0560 %modelSBML.reaction(i).notes +0561 if strfind(modelSBML.reaction(i).notes,'GENE_ASSOCIATION') +0562 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE_ASSOCIATION'); +0563 elseif strfind(modelSBML.reaction(i).notes,'GENE ASSOCIATION') +0564 geneAssociation=parseNote(modelSBML.reaction(i).notes,'GENE ASSOCIATION'); +0565 else +0566 geneAssociation=''; +0567 end +0568 if ~isempty(geneAssociation) +0569 %This adds the grRules. The gene list and rxnGeneMat are created +0570 %later +0571 grRules{counter}=geneAssociation; +0572 end +0573 end +0574 if isempty(grRules{counter}) && ~isempty(modelSBML.reaction(i).modifier) +0575 rules=''; +0576 for j=1:numel(modelSBML.reaction(i).modifier) +0577 modifier=modelSBML.reaction(i).modifier(j).species; +0578 if ~isempty(modifier) +0579 if strcmpi(modifier(1:2),'E_') +0580 index=find(strcmp(modifier,geneIDs)); +0581 %This should be unique and in the geneIDs list, +0582 %otherwise something is wrong +0583 if numel(index)~=1 +0584 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; +0585 dispEM(EM); +0586 end +0587 if ~isempty(rules) +0588 rules=[rules ' or (' geneNames{index} ')']; +0589 else +0590 rules=['(' geneNames{index} ')']; +0591 end +0592 elseif strcmp(modifier(1:2),'s_') +0593 index=find(strcmp(modifier,metaboliteIDs)); +0594 %This should be unique and in the geneIDs list, +0595 %otherwise something is wrong +0596 if numel(index)~=1 +0597 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; +0598 dispEM(EM); +0599 end +0600 if ~isempty(rules) +0601 rules=[rules ' or (' metaboliteIDs{index} ')']; +0602 else +0603 rules=['(' metaboliteIDs{index} ')']; +0604 end +0605 else +0606 %It seems to be a complex. Add the corresponding +0607 %genes from the name of the complex (not the +0608 %reaction that creates it) +0609 index=find(strcmp(modifier,complexIDs)); +0610 if numel(index)==1 +0611 if ~isempty(rules) +0612 rules=[rules ' or (' strrep(complexNames{index},':',' and ') ')']; +0613 else +0614 rules=['(' strrep(complexNames{index},':',' and ') ')']; +0615 end +0616 else +0617 %Could not find a complex +0618 EM=['Could not get the gene association data from reaction ' reactionIDs{i}]; +0619 dispEM(EM); +0620 end +0621 end +0622 end +0623 end +0624 grRules{counter}=rules; +0625 grRulesFromModifier{counter}=rules;%Backup copy for grRules, useful to parse Yeast 7.6 +0626 end +0627 +0628 %Add reaction compartment +0629 if isfield(modelSBML.reaction(i),'compartment') +0630 if ~isempty(modelSBML.reaction(i).compartment) +0631 rxnComp=modelSBML.reaction(i).compartment; +0632 else +0633 rxnComp=''; +0634 end +0635 elseif isfield(modelSBML.reaction(i),'notes') +0636 rxnComp=parseNote(modelSBML.reaction(i).notes,'COMPARTMENT'); +0637 end +0638 if ~isempty(rxnComp) +0639 %Find it in the compartment list +0640 [~, J]=ismember(rxnComp,compartmentIDs); +0641 rxnComps(counter)=J; +0642 end +0643 +0644 %Get other Miriam fields. This may include for example database indexes +0645 %to organism-specific databases. EC-codes are supported by the COBRA +0646 %Toolbox format and are therefore loaded separately +0647 if isSBML2COBRA==false +0648 miriamStruct=parseMiriam(modelSBML.reaction(i).annotation); +0649 rxnMiriams{counter}=miriamStruct; +0650 if isfield(modelSBML.reaction(i),'notes') +0651 subsystems{counter,1}=cellstr(parseNote(modelSBML.reaction(i).notes,'SUBSYSTEM')); +0652 subsystems{counter,1}(cellfun('isempty',subsystems{counter,1})) = []; +0653 if strfind(modelSBML.reaction(i).notes,'Confidence Level') +0654 rxnconfidencescores(counter)=str2num(parseNote(modelSBML.reaction(i).notes,'Confidence Level')); +0655 end +0656 rxnreferences{counter,1}=parseNote(modelSBML.reaction(i).notes,'AUTHORS'); +0657 rxnnotes{counter,1}=parseNote(modelSBML.reaction(i).notes,'NOTES'); +0658 end +0659 end +0660 +0661 %Get SBO terms +0662 if isfield(modelSBML.reaction(i),'sboTerm') +0663 rxnMiriams{counter} = addSBOtoMiriam(rxnMiriams{counter}, modelSBML.reaction(i).sboTerm); +0664 end +0665 +0666 %Get ec-codes +0667 eccode=''; +0668 if ~isempty(modelSBML.reaction(i).annotation) +0669 if strfind(modelSBML.reaction(i).annotation,'urn:miriam:ec-code') +0670 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'urn:miriam:',':','ec-code'); +0671 elseif strfind(modelSBML.reaction(i).annotation,'http://identifiers.org/ec-code') +0672 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'http://identifiers.org/','/','ec-code'); +0673 elseif strfind(modelSBML.reaction(i).annotation,'https://identifiers.org/ec-code') +0674 eccode=parseAnnotation(modelSBML.reaction(i).annotation,'https://identifiers.org/','/','ec-code'); 0675 end -0676 end -0677 eccodes{counter}=eccode; -0678 -0679 %Add all reactants -0680 for j=1:numel(modelSBML.reaction(i).reactant) -0681 %Get the index of the metabolite in metaboliteIDs. External -0682 %metabolites will be removed at a later stage -0683 metIndex=find(strcmp(modelSBML.reaction(i).reactant(j).species,metaboliteIDs),1); -0684 if isempty(metIndex) -0685 EM=['Could not find metabolite ' modelSBML.reaction(i).reactant(j).species ' in reaction ' reactionIDs{counter}]; -0686 dispEM(EM); -0687 end -0688 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).reactant(j).stoichiometry*-1; -0689 end -0690 -0691 %Add all products -0692 for j=1:numel(modelSBML.reaction(i).product) -0693 %Get the index of the metabolite in metaboliteIDs. -0694 metIndex=find(strcmp(modelSBML.reaction(i).product(j).species,metaboliteIDs),1); -0695 if isempty(metIndex) -0696 EM=['Could not find metabolite ' modelSBML.reaction(i).product(j).species ' in reaction ' reactionIDs{counter}]; -0697 dispEM(EM); -0698 end -0699 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).product(j).stoichiometry; -0700 end -0701 end -0702 -0703 %if FBC, objective function is separately defined. Multiple objective -0704 %functions can be defined, one is set as active -0705 if isfield(modelSBML, 'fbc_activeObjective') -0706 obj=modelSBML.fbc_activeObjective; -0707 for i=1:numel(modelSBML.fbc_objective) -0708 if strcmp(obj,modelSBML.fbc_objective(i).fbc_id) -0709 rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction; -0710 rxn=regexprep(rxn,'^R_',''); -0711 idx=find(ismember(reactionIDs,rxn)); -0712 reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient; -0713 end -0714 end -0715 end -0716 -0717 %subSystems can be stored as groups instead of in annotations -0718 if isfield(modelSBML,'groups_group') -0719 for i=1:numel(modelSBML.groups_group) -0720 groupreactions={modelSBML.groups_group(i).groups_member(:).groups_idRef}; -0721 groupreactions=regexprep(groupreactions,'^R_',''); -0722 [~, idx] = ismember(groupreactions, reactionIDs); -0723 if any(idx) -0724 for j=1:numel(idx) -0725 if isempty(subsystems{idx(j)}) % First subsystem -0726 subsystems{idx(j)} = {modelSBML.groups_group(i).groups_name}; -0727 else % Consecutive subsystems: concatenate -0728 subsystems{idx(j)} = horzcat(subsystems{idx(j)}, modelSBML.groups_group(i).groups_name); -0729 end -0730 end -0731 end -0732 end -0733 end -0734 -0735 %Shrink the structures if complex-forming reactions had to be skipped -0736 reactionNames=reactionNames(1:counter); -0737 reactionIDs=reactionIDs(1:counter); -0738 subsystems=subsystems(1:counter); -0739 eccodes=eccodes(1:counter); -0740 rxnconfidencescores=rxnconfidencescores(1:counter); -0741 rxnreferences=rxnreferences(1:counter); -0742 rxnnotes=rxnnotes(1:counter); -0743 grRules=grRules(1:counter); -0744 rxnMiriams=rxnMiriams(1:counter); -0745 reactionReversibility=reactionReversibility(1:counter); -0746 reactionUB=reactionUB(1:counter); -0747 reactionLB=reactionLB(1:counter); -0748 reactionObjective=reactionObjective(1:counter); -0749 S=S(:,1:counter); -0750 -0751 model.name=modelSBML.name; -0752 model.id=regexprep(modelSBML.id,'^M_',''); % COBRA adds M_ prefix -0753 model.rxns=reactionIDs; -0754 model.mets=metaboliteIDs; -0755 model.S=sparse(S); -0756 model.lb=reactionLB; -0757 model.ub=reactionUB; -0758 model.rev=reactionReversibility; -0759 model.c=reactionObjective; -0760 model.b=zeros(numel(metaboliteIDs),1); -0761 model.comps=compartmentIDs; -0762 model.compNames=compartmentNames; -0763 model.rxnConfidenceScores=rxnconfidencescores; -0764 model.rxnReferences=rxnreferences; -0765 model.rxnNotes=rxnnotes; -0766 -0767 %Load annotation if available. If there are several authors, only the first -0768 %author credentials are imported -0769 if isfield(modelSBML,'annotation') -0770 endString='</'; -0771 I=strfind(modelSBML.annotation,endString); -0772 J=strfind(modelSBML.annotation,'<vCard:Family>'); -0773 if any(J) -0774 model.annotation.familyName=modelSBML.annotation(J(1)+14:I(find(I>J(1),1))-1); -0775 end -0776 J=strfind(modelSBML.annotation,'<vCard:Given>'); -0777 if any(J) -0778 model.annotation.givenName=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); -0779 end -0780 J=strfind(modelSBML.annotation,'<vCard:EMAIL>'); -0781 if any(J) -0782 model.annotation.email=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); -0783 end -0784 J=strfind(modelSBML.annotation,'<vCard:Orgname>'); -0785 if any(J) -0786 model.annotation.organization=modelSBML.annotation(J(1)+15:I(find(I>J(1),1))-1); -0787 end -0788 endString='"/>'; -0789 I=strfind(modelSBML.annotation,endString); -0790 if strfind(modelSBML.annotation,'"urn:miriam:') -0791 J=strfind(modelSBML.annotation,'"urn:miriam:'); -0792 if any(J) -0793 model.annotation.taxonomy=modelSBML.annotation(J+12:I(find(I>J,1))-1); -0794 end -0795 else -0796 J=strfind(modelSBML.annotation,'"http://identifiers.org/'); -0797 if any(J) -0798 model.annotation.taxonomy=modelSBML.annotation(J+24:I(find(I>J,1))-1); -0799 else -0800 J=strfind(modelSBML.annotation,'"https://identifiers.org/'); -0801 if any(J) -0802 model.annotation.taxonomy=modelSBML.annotation(J+25:I(find(I>J,1))-1); -0803 end -0804 end -0805 end -0806 end -0807 if isfield(modelSBML,'notes') -0808 startString=strfind(modelSBML.notes,'xhtml">'); -0809 endString=strfind(modelSBML.notes,'</body>'); -0810 if any(startString) && any(endString) -0811 model.annotation.note=modelSBML.notes(startString+7:endString-1); -0812 model.annotation.note=regexprep(model.annotation.note,'<p>|</p>',''); -0813 model.annotation.note=strtrim(model.annotation.note); -0814 end -0815 end -0816 -0817 if any(~cellfun(@isempty,compartmentOutside)) -0818 model.compOutside=compartmentOutside; -0819 end -0820 -0821 model.rxnNames=reactionNames; -0822 model.metNames=metaboliteNames; -0823 -0824 %Match the compartments for metabolites -0825 [~, J]=ismember(metaboliteCompartments,model.comps); -0826 model.metComps=J; -0827 -0828 %If any genes have been loaded (only for the new format) -0829 if ~isempty(geneNames) -0830 %In some rare cases geneNames may not necessarily be used in grRules. -0831 %That is true for Yeast 7.6. It's therefore important to change gene -0832 %systematic names to geneIDs in sophisticated way. Gene systematic -0833 %names are not unique, since exactly the same name may be in different -0834 %compartments -0835 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) -0836 geneShortNames=geneNames; -0837 %geneShortNames contain compartments as well, so these are removed -0838 geneShortNames=regexprep(geneShortNames,' \[.+$',''); -0839 %grRules obtained from modifier fields contain geneNames. These are -0840 %changed into geneIDs. grRulesFromModifier is a good way to have -0841 %geneIDs and rxns association when it's important to resolve -0842 %systematic name ambiguities -0843 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0844 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); -0845 -0846 %Yeast 7.6 contains several metabolites, which were used in gene -0847 %associations. For that reason, the list of species ID is created -0848 %and we then check whether any of them have kegg.genes annotation -0849 %thereby obtaining systematic gene names -0850 geneShortNames=vertcat(geneShortNames,metaboliteNames); -0851 geneIDs=vertcat(geneIDs,metaboliteIDs); -0852 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); -0853 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); -0854 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); -0855 -0856 %Now we retain information for only these entries, which have -0857 %kegg.genes annotation -0858 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); -0859 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); -0860 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); -0861 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); -0862 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); -0863 %Now we reorder geneIDs and geneSystNames by geneSystNames string -0864 %length -0865 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later -0866 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); -0867 geneIDs = geneIDs(Indx); -0868 geneSystNames = geneSystNames(Indx); -0869 for i=1:numel(geneSystNames) -0870 for j=1:numel(grRules) -0871 if strfind(grRules{j},geneSystNames{i}) -0872 if ~isempty(grRules{j}) -0873 if sum(ismember(geneSystNames,geneSystNames{i}))==1 -0874 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); -0875 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 -0876 counter=0; -0877 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); -0878 for k=1:numel(ovrlpIDs) -0879 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) -0880 counter=counter+1; -0881 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); -0882 end -0883 if counter>1 -0884 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; -0885 dispEM(EM); -0886 end -0887 end -0888 end -0889 end -0890 end -0891 end -0892 end -0893 end -0894 model.genes=geneNames; -0895 model.grRules=grRules; -0896 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0897 model.grRules = grRules; -0898 model.rxnGeneMat = rxnGeneMat; -0899 -0900 %Match the compartments for genes -0901 [~, J]=ismember(geneCompartments,model.comps); -0902 model.geneComps=J; -0903 else -0904 if ~all(cellfun(@isempty,grRules)) -0905 %If fbc_geneProduct exists, follow the specified gene order, such -0906 %that matching geneShortNames in function below will work -0907 if isfield(modelSBML,'fbc_geneProduct') -0908 genes={modelSBML.fbc_geneProduct.fbc_id}; -0909 -0910 %Get gene Miriams if they were not retrieved above (this occurs -0911 %when genes are stored as fbc_geneProduct instead of species) -0912 if isempty(geneMiriams) -0913 geneMiriams = cell(numel(genes),1); -0914 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && numel(unique([modelSBML.fbc_geneProduct.sboTerm])) == 1 -0915 %If all the SBO terms are identical, don't add them to geneMiriams -0916 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); -0917 end -0918 for i = 1:numel(genes) -0919 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); -0920 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') -0921 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); -0922 end +0676 elseif isfield(modelSBML.reaction(i),'notes') +0677 if strfind(modelSBML.reaction(i).notes,'EC Number') +0678 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'EC Number')]; +0679 elseif strfind(modelSBML.reaction(i).notes,'PROTEIN_CLASS') +0680 eccode=[eccode parseNote(modelSBML.reaction(i).notes,'PROTEIN_CLASS')]; +0681 end +0682 end +0683 eccodes{counter}=eccode; +0684 +0685 %Add all reactants +0686 for j=1:numel(modelSBML.reaction(i).reactant) +0687 %Get the index of the metabolite in metaboliteIDs. External +0688 %metabolites will be removed at a later stage +0689 metIndex=find(strcmp(modelSBML.reaction(i).reactant(j).species,metaboliteIDs),1); +0690 if isempty(metIndex) +0691 EM=['Could not find metabolite ' modelSBML.reaction(i).reactant(j).species ' in reaction ' reactionIDs{counter}]; +0692 dispEM(EM); +0693 end +0694 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).reactant(j).stoichiometry*-1; +0695 end +0696 +0697 %Add all products +0698 for j=1:numel(modelSBML.reaction(i).product) +0699 %Get the index of the metabolite in metaboliteIDs. +0700 metIndex=find(strcmp(modelSBML.reaction(i).product(j).species,metaboliteIDs),1); +0701 if isempty(metIndex) +0702 EM=['Could not find metabolite ' modelSBML.reaction(i).product(j).species ' in reaction ' reactionIDs{counter}]; +0703 dispEM(EM); +0704 end +0705 S(metIndex,counter)=S(metIndex,counter)+modelSBML.reaction(i).product(j).stoichiometry; +0706 end +0707 end +0708 +0709 %if FBC, objective function is separately defined. Multiple objective +0710 %functions can be defined, one is set as active +0711 if isfield(modelSBML, 'fbc_activeObjective') +0712 obj=modelSBML.fbc_activeObjective; +0713 for i=1:numel(modelSBML.fbc_objective) +0714 if strcmp(obj,modelSBML.fbc_objective(i).fbc_id) +0715 rxn=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_reaction; +0716 rxn=regexprep(rxn,'^R_',''); +0717 idx=find(ismember(reactionIDs,rxn)); +0718 reactionObjective(idx)=modelSBML.fbc_objective(i).fbc_fluxObjective.fbc_coefficient; +0719 end +0720 end +0721 end +0722 +0723 %subSystems can be stored as groups instead of in annotations +0724 if isfield(modelSBML,'groups_group') +0725 for i=1:numel(modelSBML.groups_group) +0726 groupreactions={modelSBML.groups_group(i).groups_member(:).groups_idRef}; +0727 groupreactions=regexprep(groupreactions,'^R_',''); +0728 [~, idx] = ismember(groupreactions, reactionIDs); +0729 if any(idx) +0730 for j=1:numel(idx) +0731 if isempty(subsystems{idx(j)}) % First subsystem +0732 subsystems{idx(j)} = {modelSBML.groups_group(i).groups_name}; +0733 else % Consecutive subsystems: concatenate +0734 subsystems{idx(j)} = horzcat(subsystems{idx(j)}, modelSBML.groups_group(i).groups_name); +0735 end +0736 end +0737 end +0738 end +0739 end +0740 +0741 %Shrink the structures if complex-forming reactions had to be skipped +0742 reactionNames=reactionNames(1:counter); +0743 reactionIDs=reactionIDs(1:counter); +0744 subsystems=subsystems(1:counter); +0745 eccodes=eccodes(1:counter); +0746 rxnconfidencescores=rxnconfidencescores(1:counter); +0747 rxnreferences=rxnreferences(1:counter); +0748 rxnnotes=rxnnotes(1:counter); +0749 grRules=grRules(1:counter); +0750 rxnMiriams=rxnMiriams(1:counter); +0751 reactionReversibility=reactionReversibility(1:counter); +0752 reactionUB=reactionUB(1:counter); +0753 reactionLB=reactionLB(1:counter); +0754 reactionObjective=reactionObjective(1:counter); +0755 S=S(:,1:counter); +0756 +0757 model.name=modelSBML.name; +0758 model.id=regexprep(modelSBML.id,'^M_',''); % COBRA adds M_ prefix +0759 model.rxns=reactionIDs; +0760 model.mets=metaboliteIDs; +0761 model.S=sparse(S); +0762 model.lb=reactionLB; +0763 model.ub=reactionUB; +0764 model.rev=reactionReversibility; +0765 model.c=reactionObjective; +0766 model.b=zeros(numel(metaboliteIDs),1); +0767 model.comps=compartmentIDs; +0768 model.compNames=compartmentNames; +0769 model.rxnConfidenceScores=rxnconfidencescores; +0770 model.rxnReferences=rxnreferences; +0771 model.rxnNotes=rxnnotes; +0772 +0773 %Load annotation if available. If there are several authors, only the first +0774 %author credentials are imported +0775 if isfield(modelSBML,'annotation') +0776 endString='</'; +0777 I=strfind(modelSBML.annotation,endString); +0778 J=strfind(modelSBML.annotation,'<vCard:Family>'); +0779 if any(J) +0780 model.annotation.familyName=modelSBML.annotation(J(1)+14:I(find(I>J(1),1))-1); +0781 end +0782 J=strfind(modelSBML.annotation,'<vCard:Given>'); +0783 if any(J) +0784 model.annotation.givenName=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); +0785 end +0786 J=strfind(modelSBML.annotation,'<vCard:EMAIL>'); +0787 if any(J) +0788 model.annotation.email=modelSBML.annotation(J(1)+13:I(find(I>J(1),1))-1); +0789 end +0790 J=strfind(modelSBML.annotation,'<vCard:Orgname>'); +0791 if any(J) +0792 model.annotation.organization=modelSBML.annotation(J(1)+15:I(find(I>J(1),1))-1); +0793 end +0794 endString='"/>'; +0795 I=strfind(modelSBML.annotation,endString); +0796 if strfind(modelSBML.annotation,'"urn:miriam:') +0797 J=strfind(modelSBML.annotation,'"urn:miriam:'); +0798 if any(J) +0799 model.annotation.taxonomy=modelSBML.annotation(J+12:I(find(I>J,1))-1); +0800 end +0801 else +0802 J=strfind(modelSBML.annotation,'"http://identifiers.org/'); +0803 if any(J) +0804 model.annotation.taxonomy=modelSBML.annotation(J+24:I(find(I>J,1))-1); +0805 else +0806 J=strfind(modelSBML.annotation,'"https://identifiers.org/'); +0807 if any(J) +0808 model.annotation.taxonomy=modelSBML.annotation(J+25:I(find(I>J,1))-1); +0809 end +0810 end +0811 end +0812 end +0813 if isfield(modelSBML,'notes') +0814 startString=strfind(modelSBML.notes,'xhtml">'); +0815 endString=strfind(modelSBML.notes,'</body>'); +0816 if any(startString) && any(endString) +0817 model.annotation.note=modelSBML.notes(startString+7:endString-1); +0818 model.annotation.note=regexprep(model.annotation.note,'<p>|</p>',''); +0819 model.annotation.note=strtrim(model.annotation.note); +0820 end +0821 end +0822 +0823 if any(~cellfun(@isempty,compartmentOutside)) +0824 model.compOutside=compartmentOutside; +0825 end +0826 +0827 model.rxnNames=reactionNames; +0828 model.metNames=metaboliteNames; +0829 +0830 %Match the compartments for metabolites +0831 [~, J]=ismember(metaboliteCompartments,model.comps); +0832 model.metComps=J; +0833 +0834 %If any genes have been loaded (only for the new format) +0835 if ~isempty(geneNames) +0836 %In some rare cases geneNames may not necessarily be used in grRules. +0837 %That is true for Yeast 7.6. It's therefore important to change gene +0838 %systematic names to geneIDs in sophisticated way. Gene systematic +0839 %names are not unique, since exactly the same name may be in different +0840 %compartments +0841 if all(cellfun(@isempty,strfind(grRules,geneNames{1}))) +0842 geneShortNames=geneNames; +0843 %geneShortNames contain compartments as well, so these are removed +0844 geneShortNames=regexprep(geneShortNames,' \[.+$',''); +0845 %grRules obtained from modifier fields contain geneNames. These are +0846 %changed into geneIDs. grRulesFromModifier is a good way to have +0847 %geneIDs and rxns association when it's important to resolve +0848 %systematic name ambiguities +0849 grRulesFromModifier=regexprep(regexprep(grRulesFromModifier,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0850 grRules=regexprep(regexprep(grRules,'\[|\]','_'),regexprep(geneNames,'\[|\]','_'),geneIDs); +0851 +0852 %Yeast 7.6 contains several metabolites, which were used in gene +0853 %associations. For that reason, the list of species ID is created +0854 %and we then check whether any of them have kegg.genes annotation +0855 %thereby obtaining systematic gene names +0856 geneShortNames=vertcat(geneShortNames,metaboliteNames); +0857 geneIDs=vertcat(geneIDs,metaboliteIDs); +0858 geneSystNames=extractMiriam(vertcat(geneMiriams,metaboliteMiriams),'kegg.genes'); +0859 geneCompartments=vertcat(geneCompartments,metaboliteCompartments); +0860 geneMiriams=vertcat(geneMiriams,metaboliteMiriams); +0861 +0862 %Now we retain information for only these entries, which have +0863 %kegg.genes annotation +0864 geneShortNames=geneShortNames(~cellfun('isempty',geneSystNames)); +0865 geneIDs=geneIDs(~cellfun('isempty',geneSystNames)); +0866 geneSystNames=geneSystNames(~cellfun('isempty',geneSystNames)); +0867 geneCompartments=geneCompartments(~cellfun('isempty',geneSystNames)); +0868 geneMiriams=geneMiriams(~cellfun('isempty',geneSystNames)); +0869 %Now we reorder geneIDs and geneSystNames by geneSystNames string +0870 %length +0871 geneNames=geneIDs;%Backuping geneIDs, since we need unsorted order for later +0872 [~, Indx] = sort(cellfun('size', geneSystNames, 2), 'descend'); +0873 geneIDs = geneIDs(Indx); +0874 geneSystNames = geneSystNames(Indx); +0875 for i=1:numel(geneSystNames) +0876 for j=1:numel(grRules) +0877 if strfind(grRules{j},geneSystNames{i}) +0878 if ~isempty(grRules{j}) +0879 if sum(ismember(geneSystNames,geneSystNames{i}))==1 +0880 grRules{j}=regexprep(grRules{j},geneSystNames{i},geneIDs{i}); +0881 elseif sum(ismember(geneSystNames,geneSystNames{i}))>1 +0882 counter=0; +0883 ovrlpIDs=geneIDs(ismember(geneSystNames,geneSystNames{i})); +0884 for k=1:numel(ovrlpIDs) +0885 if strfind(grRulesFromModifier{j},ovrlpIDs{k}) +0886 counter=counter+1; +0887 grRules{j}=regexprep(grRules{j},geneSystNames{i},ovrlpIDs{k}); +0888 end +0889 if counter>1 +0890 EM=['Gene association is ambiguous for reaction ' modelSBML.reaction(j).id]; +0891 dispEM(EM); +0892 end +0893 end +0894 end +0895 end +0896 end +0897 end +0898 end +0899 end +0900 model.genes=geneNames; +0901 model.grRules=grRules; +0902 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0903 model.grRules = grRules; +0904 model.rxnGeneMat = rxnGeneMat; +0905 +0906 %Match the compartments for genes +0907 [~, J]=ismember(geneCompartments,model.comps); +0908 model.geneComps=J; +0909 else +0910 if ~all(cellfun(@isempty,grRules)) +0911 %If fbc_geneProduct exists, follow the specified gene order, such +0912 %that matching geneShortNames in function below will work +0913 if isfield(modelSBML,'fbc_geneProduct') +0914 genes={modelSBML.fbc_geneProduct.fbc_id}; +0915 +0916 %Get gene Miriams if they were not retrieved above (this occurs +0917 %when genes are stored as fbc_geneProduct instead of species) +0918 if isempty(geneMiriams) +0919 geneMiriams = cell(numel(genes),1); +0920 if isfield(modelSBML.fbc_geneProduct,'sboTerm') && numel(unique([modelSBML.fbc_geneProduct.sboTerm])) == 1 +0921 %If all the SBO terms are identical, don't add them to geneMiriams +0922 modelSBML.fbc_geneProduct = rmfield(modelSBML.fbc_geneProduct,'sboTerm'); 0923 end -0924 end -0925 else -0926 genes=getGeneList(grRules); -0927 end -0928 if strcmpi(genes{1}(1:2),'G_') -0929 genes=regexprep(genes,'^G_',''); -0930 grRules=regexprep(grRules,'^G_',''); -0931 grRules=regexprep(grRules,'\(G_','('); -0932 grRules=regexprep(grRules,' G_',' '); +0924 for i = 1:numel(genes) +0925 geneMiriams{i}=parseMiriam(modelSBML.fbc_geneProduct(i).annotation); +0926 if isfield(modelSBML.fbc_geneProduct(i),'sboTerm') +0927 geneMiriams{i} = addSBOtoMiriam(geneMiriams{i},modelSBML.fbc_geneProduct(i).sboTerm); +0928 end +0929 end +0930 end +0931 else +0932 genes=getGeneList(grRules); 0933 end -0934 model.genes=genes; -0935 model.grRules=grRules; -0936 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0937 model.grRules = grRules; -0938 model.rxnGeneMat = rxnGeneMat; -0939 end -0940 end -0941 -0942 if all(cellfun(@isempty,geneShortNames)) -0943 if isfield(modelSBML,'fbc_geneProduct') -0944 for i=1:numel(genes) -0945 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) -0946 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; -0947 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) -0948 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; -0949 else -0950 geneShortNames{i,1}=''; -0951 end -0952 end -0953 end -0954 end -0955 -0956 %If any InChIs have been loaded -0957 if any(~cellfun(@isempty,metaboliteInChI)) -0958 model.inchis=metaboliteInChI; -0959 end -0960 -0961 %If any formulas have been loaded -0962 if any(~cellfun(@isempty,metaboliteFormula)) -0963 model.metFormulas=metaboliteFormula; -0964 end -0965 -0966 %If any charges have been loaded -0967 if ~isempty(metaboliteCharges) -0968 model.metCharges=metaboliteCharges; -0969 end -0970 -0971 %If any gene short names have been loaded -0972 if any(~cellfun(@isempty,geneShortNames)) -0973 model.geneShortNames=geneShortNames; -0974 end -0975 -0976 %If any Miriam strings for compartments have been loaded -0977 if any(~cellfun(@isempty,compartmentMiriams)) -0978 model.compMiriams=compartmentMiriams; -0979 end -0980 -0981 %If any Miriam strings for metabolites have been loaded -0982 if any(~cellfun(@isempty,metaboliteMiriams)) -0983 model.metMiriams=metaboliteMiriams; -0984 end -0985 -0986 %If any subsystems have been loaded -0987 if any(~cellfun(@isempty,subsystems)) -0988 model.subSystems=subsystems; -0989 end -0990 if any(rxnComps) -0991 if all(rxnComps) -0992 model.rxnComps=rxnComps; -0993 else -0994 if supressWarnings==false -0995 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; -0996 dispEM(EM,false,model.rxns(rxnComps==0)); -0997 end -0998 end -0999 end -1000 -1001 %If any ec-codes have been loaded -1002 if any(~cellfun(@isempty,eccodes)) -1003 model.eccodes=eccodes; -1004 end -1005 -1006 %If any Miriam strings for reactions have been loaded -1007 if any(~cellfun(@isempty,rxnMiriams)) -1008 model.rxnMiriams=rxnMiriams; -1009 end -1010 -1011 %If any Miriam strings for genes have been loaded -1012 if any(~cellfun(@isempty,geneMiriams)) -1013 model.geneMiriams=geneMiriams; -1014 end -1015 -1016 model.unconstrained=metaboliteUnconstrained; -1017 -1018 %Convert SBML IDs back into their original strings. Here we are using part -1019 %from convertSBMLID, originating from the COBRA Toolbox -1020 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); -1021 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); -1022 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); -1023 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); -1024 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); -1025 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); -1026 -1027 %Remove unused fields -1028 if isempty(model.annotation) -1029 model=rmfield(model,'annotation'); -1030 end -1031 if isempty(model.compOutside) -1032 model=rmfield(model,'compOutside'); -1033 end -1034 if isempty(model.compMiriams) -1035 model=rmfield(model,'compMiriams'); +0934 if strcmpi(genes{1}(1:2),'G_') +0935 genes=regexprep(genes,'^G_',''); +0936 grRules=regexprep(grRules,'^G_',''); +0937 grRules=regexprep(grRules,'\(G_','('); +0938 grRules=regexprep(grRules,' G_',' '); +0939 end +0940 model.genes=genes; +0941 model.grRules=grRules; +0942 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0943 model.grRules = grRules; +0944 model.rxnGeneMat = rxnGeneMat; +0945 end +0946 end +0947 +0948 if all(cellfun(@isempty,geneShortNames)) +0949 if isfield(modelSBML,'fbc_geneProduct') +0950 for i=1:numel(genes) +0951 if ~isempty(modelSBML.fbc_geneProduct(i).fbc_label) +0952 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_label; +0953 elseif ~isempty(modelSBML.fbc_geneProduct(i).fbc_name) +0954 geneShortNames{i,1}=modelSBML.fbc_geneProduct(i).fbc_name; +0955 else +0956 geneShortNames{i,1}=''; +0957 end +0958 end +0959 end +0960 end +0961 +0962 %If any InChIs have been loaded +0963 if any(~cellfun(@isempty,metaboliteInChI)) +0964 model.inchis=metaboliteInChI; +0965 end +0966 +0967 %If any formulas have been loaded +0968 if any(~cellfun(@isempty,metaboliteFormula)) +0969 model.metFormulas=metaboliteFormula; +0970 end +0971 +0972 %If any charges have been loaded +0973 if ~isempty(metaboliteCharges) +0974 model.metCharges=metaboliteCharges; +0975 end +0976 +0977 %If any gene short names have been loaded +0978 if any(~cellfun(@isempty,geneShortNames)) +0979 model.geneShortNames=geneShortNames; +0980 end +0981 +0982 %If any Miriam strings for compartments have been loaded +0983 if any(~cellfun(@isempty,compartmentMiriams)) +0984 model.compMiriams=compartmentMiriams; +0985 end +0986 +0987 %If any Miriam strings for metabolites have been loaded +0988 if any(~cellfun(@isempty,metaboliteMiriams)) +0989 model.metMiriams=metaboliteMiriams; +0990 end +0991 +0992 %If any subsystems have been loaded +0993 if any(~cellfun(@isempty,subsystems)) +0994 model.subSystems=subsystems; +0995 end +0996 if any(rxnComps) +0997 if all(rxnComps) +0998 model.rxnComps=rxnComps; +0999 else +1000 if supressWarnings==false +1001 EM='The compartments for the following reactions could not be matched. Ignoring reaction compartment information'; +1002 dispEM(EM,false,model.rxns(rxnComps==0)); +1003 end +1004 end +1005 end +1006 +1007 %If any ec-codes have been loaded +1008 if any(~cellfun(@isempty,eccodes)) +1009 model.eccodes=eccodes; +1010 end +1011 +1012 %If any Miriam strings for reactions have been loaded +1013 if any(~cellfun(@isempty,rxnMiriams)) +1014 model.rxnMiriams=rxnMiriams; +1015 end +1016 +1017 %If any Miriam strings for genes have been loaded +1018 if any(~cellfun(@isempty,geneMiriams)) +1019 model.geneMiriams=geneMiriams; +1020 end +1021 +1022 model.unconstrained=metaboliteUnconstrained; +1023 +1024 %Convert SBML IDs back into their original strings. Here we are using part +1025 %from convertSBMLID, originating from the COBRA Toolbox +1026 model.rxns=regexprep(model.rxns,'__([0-9]+)__','${char(str2num($1))}'); +1027 model.mets=regexprep(model.mets,'__([0-9]+)__','${char(str2num($1))}'); +1028 model.comps=regexprep(model.comps,'__([0-9]+)__','${char(str2num($1))}'); +1029 model.grRules=regexprep(model.grRules,'__([0-9]+)__','${char(str2num($1))}'); +1030 model.genes=regexprep(model.genes,'__([0-9]+)__','${char(str2num($1))}'); +1031 model.id=regexprep(model.id,'__([0-9]+)__','${char(str2num($1))}'); +1032 +1033 %Remove unused fields +1034 if isempty(model.annotation) +1035 model=rmfield(model,'annotation'); 1036 end -1037 if isempty(model.rxnComps) -1038 model=rmfield(model,'rxnComps'); +1037 if isempty(model.compOutside) +1038 model=rmfield(model,'compOutside'); 1039 end -1040 if isempty(model.grRules) -1041 model=rmfield(model,'grRules'); +1040 if isempty(model.compMiriams) +1041 model=rmfield(model,'compMiriams'); 1042 end -1043 if isempty(model.rxnGeneMat) -1044 model=rmfield(model,'rxnGeneMat'); +1043 if isempty(model.rxnComps) +1044 model=rmfield(model,'rxnComps'); 1045 end -1046 if isempty(model.subSystems) -1047 model=rmfield(model,'subSystems'); -1048 else -1049 model.subSystems(cellfun(@isempty,subsystems))={{''}}; -1050 end -1051 if isempty(model.eccodes) -1052 model=rmfield(model,'eccodes'); -1053 end -1054 if isempty(model.rxnMiriams) -1055 model=rmfield(model,'rxnMiriams'); +1046 if isempty(model.grRules) +1047 model=rmfield(model,'grRules'); +1048 end +1049 if isempty(model.rxnGeneMat) +1050 model=rmfield(model,'rxnGeneMat'); +1051 end +1052 if isempty(model.subSystems) +1053 model=rmfield(model,'subSystems'); +1054 else +1055 model.subSystems(cellfun(@isempty,subsystems))={{''}}; 1056 end -1057 if cellfun(@isempty,model.rxnNotes) -1058 model=rmfield(model,'rxnNotes'); +1057 if isempty(model.eccodes) +1058 model=rmfield(model,'eccodes'); 1059 end -1060 if cellfun(@isempty,model.rxnReferences) -1061 model=rmfield(model,'rxnReferences'); +1060 if isempty(model.rxnMiriams) +1061 model=rmfield(model,'rxnMiriams'); 1062 end -1063 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) -1064 model=rmfield(model,'rxnConfidenceScores'); +1063 if cellfun(@isempty,model.rxnNotes) +1064 model=rmfield(model,'rxnNotes'); 1065 end -1066 if isempty(model.genes) -1067 model=rmfield(model,'genes'); -1068 elseif isrow(model.genes) -1069 model.genes=transpose(model.genes); -1070 end -1071 if isempty(model.geneComps) -1072 model=rmfield(model,'geneComps'); -1073 end -1074 if isempty(model.geneMiriams) -1075 model=rmfield(model,'geneMiriams'); +1066 if cellfun(@isempty,model.rxnReferences) +1067 model=rmfield(model,'rxnReferences'); +1068 end +1069 if isempty(model.rxnConfidenceScores) || all(isnan(model.rxnConfidenceScores)) +1070 model=rmfield(model,'rxnConfidenceScores'); +1071 end +1072 if isempty(model.genes) +1073 model=rmfield(model,'genes'); +1074 elseif isrow(model.genes) +1075 model.genes=transpose(model.genes); 1076 end -1077 if isempty(model.geneShortNames) -1078 model=rmfield(model,'geneShortNames'); +1077 if isempty(model.geneComps) +1078 model=rmfield(model,'geneComps'); 1079 end -1080 if isempty(model.inchis) -1081 model=rmfield(model,'inchis'); +1080 if isempty(model.geneMiriams) +1081 model=rmfield(model,'geneMiriams'); 1082 end -1083 if isempty(model.metFormulas) -1084 model=rmfield(model,'metFormulas'); +1083 if isempty(model.geneShortNames) +1084 model=rmfield(model,'geneShortNames'); 1085 end -1086 if isempty(model.metMiriams) -1087 model=rmfield(model,'metMiriams'); +1086 if isempty(model.inchis) +1087 model=rmfield(model,'inchis'); 1088 end -1089 if ~any(model.metCharges) -1090 model=rmfield(model,'metCharges'); +1089 if isempty(model.metFormulas) +1090 model=rmfield(model,'metFormulas'); 1091 end -1092 -1093 %This just removes the grRules if no genes have been loaded -1094 if ~isfield(model,'genes') && isfield(model,'grRules') -1095 model=rmfield(model,'grRules'); -1096 end -1097 -1098 %Print warnings about bad structure -1099 if supressWarnings==false -1100 checkModelStruct(model,false); -1101 end -1102 -1103 if removeExcMets==true -1104 model=simplifyModel(model); -1105 end -1106 end -1107 -1108 function matchGenes=getGeneList(grRules) -1109 %Constructs the list of unique genes from grRules -1110 -1111 %Assumes that everything that isn't a paranthesis, " AND " or " or " is a -1112 %gene name -1113 genes=strrep(grRules,'(',''); -1114 genes=strrep(genes,')',''); -1115 genes=strrep(genes,' or ',' '); -1116 genes=strrep(genes,' and ',' '); -1117 genes=strrep(genes,' OR ',' '); -1118 genes=strrep(genes,' AND ',' '); -1119 genes=regexp(genes,' ','split'); -1120 -1121 allNames={}; -1122 for i=1:numel(genes) -1123 allNames=[allNames genes{i}]; -1124 end -1125 matchGenes=unique(allNames)'; +1092 if isempty(model.metMiriams) +1093 model=rmfield(model,'metMiriams'); +1094 end +1095 if ~any(model.metCharges) +1096 model=rmfield(model,'metCharges'); +1097 end +1098 +1099 %This just removes the grRules if no genes have been loaded +1100 if ~isfield(model,'genes') && isfield(model,'grRules') +1101 model=rmfield(model,'grRules'); +1102 end +1103 +1104 %Print warnings about bad structure +1105 if supressWarnings==false +1106 checkModelStruct(model,false); +1107 end +1108 +1109 if removeExcMets==true +1110 model=simplifyModel(model); +1111 end +1112 end +1113 +1114 function matchGenes=getGeneList(grRules) +1115 %Constructs the list of unique genes from grRules +1116 +1117 %Assumes that everything that isn't a paranthesis, " AND " or " or " is a +1118 %gene name +1119 genes=strrep(grRules,'(',''); +1120 genes=strrep(genes,')',''); +1121 genes=strrep(genes,' or ',' '); +1122 genes=strrep(genes,' and ',' '); +1123 genes=strrep(genes,' OR ',' '); +1124 genes=strrep(genes,' AND ',' '); +1125 genes=regexp(genes,' ','split'); 1126 -1127 %Remove the empty element if present -1128 if isempty(matchGenes{1}) -1129 matchGenes(1)=[]; +1127 allNames={}; +1128 for i=1:numel(genes) +1129 allNames=[allNames genes{i}]; 1130 end -1131 end +1131 matchGenes=unique(allNames)'; 1132 -1133 function fieldContent=parseNote(searchString,fieldName) -1134 %The function obtains the particular information from 'notes' field, using -1135 %fieldName as the dummy string -1136 -1137 fieldContent=''; +1133 %Remove the empty element if present +1134 if isempty(matchGenes{1}) +1135 matchGenes(1)=[]; +1136 end +1137 end 1138 -1139 if strfind(searchString,fieldName) -1140 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); -1141 targetString=regexprep(targetString,'<p>|</p>',''); -1142 targetString=regexprep(targetString,[fieldName, ':'],''); -1143 for i=1:numel(targetString) -1144 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1145 end -1146 fieldContent=regexprep(fieldContent,'^;|;$',''); -1147 else -1148 fieldContent=''; -1149 end -1150 end -1151 -1152 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) -1153 -1154 fieldContent=''; -1155 -1156 %Removing whitespace characters from the ending strings, which may occur in -1157 %several cases -1158 searchString=regexprep(searchString,'" />','"/>'); -1159 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); -1160 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1161 targetString=regexprep(targetString,startString,''); -1162 targetString=regexprep(targetString,[fieldName midString],''); -1163 -1164 for i=1:numel(targetString) -1165 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; -1166 end -1167 -1168 fieldContent=regexprep(fieldContent,'^;|;$',''); -1169 end -1170 -1171 function miriamStruct=parseMiriam(searchString) -1172 %Generates miriam structure from annotation field +1139 function fieldContent=parseNote(searchString,fieldName) +1140 %The function obtains the particular information from 'notes' field, using +1141 %fieldName as the dummy string +1142 +1143 fieldContent=''; +1144 +1145 if strfind(searchString,fieldName) +1146 [~,targetString] = regexp(searchString,['<p>' fieldName '.*?</p>'],'tokens','match'); +1147 targetString=regexprep(targetString,'<p>|</p>',''); +1148 targetString=regexprep(targetString,[fieldName, ':'],''); +1149 for i=1:numel(targetString) +1150 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1151 end +1152 fieldContent=regexprep(fieldContent,'^;|;$',''); +1153 else +1154 fieldContent=''; +1155 end +1156 end +1157 +1158 function fieldContent=parseAnnotation(searchString,startString,midString,fieldName) +1159 +1160 fieldContent=''; +1161 +1162 %Removing whitespace characters from the ending strings, which may occur in +1163 %several cases +1164 searchString=regexprep(searchString,'" />','"/>'); +1165 [~,targetString] = regexp(searchString,['<rdf:li rdf:resource="' startString fieldName midString '.*?"/>'],'tokens','match'); +1166 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1167 targetString=regexprep(targetString,startString,''); +1168 targetString=regexprep(targetString,[fieldName midString],''); +1169 +1170 for i=1:numel(targetString) +1171 fieldContent=[fieldContent ';' strtrim(targetString{1,i})]; +1172 end 1173 -1174 %Finding whether miriams are written in the old or the new way -1175 if strfind(searchString,'urn:miriam:') -1176 startString='urn:miriam:'; -1177 midString=':'; -1178 elseif strfind(searchString,'http://identifiers.org/') -1179 startString='http://identifiers.org/'; -1180 midString='/'; -1181 elseif strfind(searchString,'https://identifiers.org/') -1182 startString='https://identifiers.org/'; -1183 midString='/'; -1184 else -1185 miriamStruct=[]; -1186 return; -1187 end -1188 -1189 miriamStruct=[]; -1190 -1191 searchString=regexprep(searchString,'" />','"/>'); -1192 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); -1193 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); -1194 targetString=regexprep(targetString,startString,''); -1195 targetString=regexprep(targetString,midString,'/','once'); +1174 fieldContent=regexprep(fieldContent,'^;|;$',''); +1175 end +1176 +1177 function miriamStruct=parseMiriam(searchString) +1178 %Generates miriam structure from annotation field +1179 +1180 %Finding whether miriams are written in the old or the new way +1181 if strfind(searchString,'urn:miriam:') +1182 startString='urn:miriam:'; +1183 midString=':'; +1184 elseif strfind(searchString,'http://identifiers.org/') +1185 startString='http://identifiers.org/'; +1186 midString='/'; +1187 elseif strfind(searchString,'https://identifiers.org/') +1188 startString='https://identifiers.org/'; +1189 midString='/'; +1190 else +1191 miriamStruct=[]; +1192 return; +1193 end +1194 +1195 miriamStruct=[]; 1196 -1197 counter=0; -1198 for i=1:numel(targetString) -1199 if isempty(regexp(targetString{1,i},'inchi|ec-code|sbo', 'once')) -1200 counter=counter+1; -1201 miriamStruct.name{counter,1} = regexprep(targetString{1,i},'/.+','','once'); -1202 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} '/'],'','once'); -1203 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); -1204 end -1205 end -1206 end -1207 -1208 function miriam = addSBOtoMiriam(miriam,sboTerm) -1209 %Appends SBO term to miriam structure -1210 -1211 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format -1212 if isempty(miriam) -1213 miriam.name = {'sbo'}; -1214 miriam.value = sboTerm; -1215 else -1216 miriam.name(end+1) = {'sbo'}; -1217 miriam.value(end+1) = sboTerm; -1218 end -1219 -1220 end +1197 searchString=regexprep(searchString,'" />','"/>'); +1198 [~,targetString] = regexp(searchString,'<rdf:li rdf:resource=".*?"/>','tokens','match'); +1199 targetString=regexprep(targetString,'<rdf:li rdf:resource="|"/>',''); +1200 targetString=regexprep(targetString,startString,''); +1201 targetString=regexprep(targetString,midString,'/','once'); +1202 +1203 counter=0; +1204 for i=1:numel(targetString) +1205 if isempty(regexp(targetString{1,i},'inchi|ec-code|sbo', 'once')) +1206 counter=counter+1; +1207 miriamStruct.name{counter,1} = regexprep(targetString{1,i},'/.+','','once'); +1208 miriamStruct.value{counter,1} = regexprep(targetString{1,i},[miriamStruct.name{counter,1} '/'],'','once'); +1209 miriamStruct.name{counter,1} = regexprep(miriamStruct.name{counter,1},'^obo\.',''); +1210 end +1211 end +1212 end +1213 +1214 function miriam = addSBOtoMiriam(miriam,sboTerm) +1215 %Appends SBO term to miriam structure +1216 +1217 sboTerm = {['SBO:' sprintf('%07u',sboTerm)]}; % convert to proper format +1218 if isempty(miriam) +1219 miriam.name = {'sbo'}; +1220 miriam.value = sboTerm; +1221 else +1222 miriam.name(end+1) = {'sbo'}; +1223 miriam.value(end+1) = sboTerm; +1224 end +1225 +1226 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/readYAMLmodel.html b/doc/io/readYAMLmodel.html index a7d7c9c4..2914b6d9 100644 --- a/doc/io/readYAMLmodel.html +++ b/doc/io/readYAMLmodel.html @@ -635,93 +635,94 @@

    SOURCE CODE ^% end 0579 % end 0580 -0581 % Make rxnGeneMat fields -0582 [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules, model.genes); -0583 if isequal(sort(genes), sort(model.genes)) -0584 model.rxnGeneMat = rxnGeneMat; -0585 model.genes = genes; -0586 else -0587 error('The gene list and grRules are inconsistent.'); +0581 % Make rxnGeneMat fields and map to the existing model.genes field +0582 [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules); +0583 model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes)); +0584 [~,geneOrder] = ismember(genes,model.genes); +0585 if any(geneOrder == 0) +0586 error(['The grRules includes the following gene(s), that are not in '... +0587 'the list of model genes: ', genes{~geneOrder}]) 0588 end -0589 -0590 % Finalize GECKO model -0591 if isGECKO -0592 % Fill in empty fields and empty entries -0593 for i={'kcat','source','notes','eccodes'} % Even keep empty -0594 model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true); -0595 end -0596 for i={'enzymes','mw','sequence'} -0597 model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true); -0598 end -0599 model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true); -0600 model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true); -0601 % Change string to double -0602 for i={'kcat','mw','concs'} -0603 if isfield(model.ec,i{1}) -0604 model.ec.(i{1}) = str2double(model.ec.(i{1})); -0605 end -0606 end -0607 % Fill rxnEnzMat -0608 rxnIdx = cellfun(@isempty, enzStoich(:,1)); -0609 enzStoich(rxnIdx,:) = ''; -0610 rxnIdx = cell2mat(enzStoich(:,1)); -0611 [~,enzIdx] = ismember(enzStoich(:,2),model.ec.enzymes); -0612 coeffs = cell2mat(enzStoich(:,3)); -0613 model.ec.rxnEnzMat = zeros(max(rxnIdx), max(enzIdx)); -0614 linearIndices = sub2ind([max(rxnIdx), max(enzIdx)], rxnIdx, enzIdx); -0615 model.ec.rxnEnzMat(linearIndices) = coeffs; -0616 %Parse ec-codes -0617 if ~isempty(ecGecko) -0618 locs = cell2mat(ecGecko(:,1)); -0619 for i=unique(locs)' -0620 ecGeckoCat=strjoin(ecGecko(locs==i,2),';'); -0621 model.ec.eccodes{i,1}=ecGeckoCat; -0622 end -0623 emptyEc=cellfun(@isempty,model.ec.eccodes); -0624 model.ec.eccodes(emptyEc)={''}; -0625 end -0626 end -0627 -0628 if verbose -0629 fprintf(' Done!\n'); -0630 end +0589 model.rxnGeneMat(:,geneOrder) = rxnGeneMat; +0590 +0591 % Finalize GECKO model +0592 if isGECKO +0593 % Fill in empty fields and empty entries +0594 for i={'kcat','source','notes','eccodes'} % Even keep empty +0595 model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true); +0596 end +0597 for i={'enzymes','mw','sequence'} +0598 model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true); +0599 end +0600 model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true); +0601 model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true); +0602 % Change string to double +0603 for i={'kcat','mw','concs'} +0604 if isfield(model.ec,i{1}) +0605 model.ec.(i{1}) = str2double(model.ec.(i{1})); +0606 end +0607 end +0608 % Fill rxnEnzMat +0609 rxnIdx = cellfun(@isempty, enzStoich(:,1)); +0610 enzStoich(rxnIdx,:) = ''; +0611 rxnIdx = cell2mat(enzStoich(:,1)); +0612 [~,enzIdx] = ismember(enzStoich(:,2),model.ec.enzymes); +0613 coeffs = cell2mat(enzStoich(:,3)); +0614 model.ec.rxnEnzMat = zeros(max(rxnIdx), max(enzIdx)); +0615 linearIndices = sub2ind([max(rxnIdx), max(enzIdx)], rxnIdx, enzIdx); +0616 model.ec.rxnEnzMat(linearIndices) = coeffs; +0617 %Parse ec-codes +0618 if ~isempty(ecGecko) +0619 locs = cell2mat(ecGecko(:,1)); +0620 for i=unique(locs)' +0621 ecGeckoCat=strjoin(ecGecko(locs==i,2),';'); +0622 model.ec.eccodes{i,1}=ecGeckoCat; +0623 end +0624 emptyEc=cellfun(@isempty,model.ec.eccodes); +0625 model.ec.eccodes(emptyEc)={''}; +0626 end +0627 end +0628 +0629 if verbose +0630 fprintf(' Done!\n'); 0631 end -0632 -0633 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty) -0634 if nargin<5 -0635 keepEmpty=false; -0636 end -0637 if isnumeric(emptyEntry) -0638 emptyCells=isempty(model.(field)); -0639 else -0640 emptyCells=cellfun(@isempty,model.(field)); -0641 end -0642 if all(emptyCells) && ~keepEmpty -0643 model = rmfield(model, field); -0644 elseif numel(model.(field))<numel(model.(type)) -0645 model.(field)(end+1:numel(model.(type)),1)=emptyEntry; -0646 end +0632 end +0633 +0634 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty) +0635 if nargin<5 +0636 keepEmpty=false; +0637 end +0638 if isnumeric(emptyEntry) +0639 emptyCells=isempty(model.(field)); +0640 else +0641 emptyCells=cellfun(@isempty,model.(field)); +0642 end +0643 if all(emptyCells) && ~keepEmpty +0644 model = rmfield(model, field); +0645 elseif numel(model.(field))<numel(model.(type)) +0646 model.(field)(end+1:numel(model.(type)),1)=emptyEntry; 0647 end -0648 -0649 function model = readFieldValue(model, fieldName, value, pos) -0650 if numel(model.(fieldName))<pos-1 -0651 model.(fieldName)(end+1:pos,1) = {''}; -0652 end -0653 model.(fieldName)(pos,1) = {value}; -0654 end -0655 -0656 function [miriams, miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber) -0657 if isempty(key) -0658 key=miriamKey; -0659 else -0660 miriamKey=key; -0661 end -0662 if ~isempty(value) -0663 miriams(entryNumber,1:3) = {pos, key, strip(value)}; -0664 else -0665 entryNumber = entryNumber - 1; -0666 end -0667 end +0648 end +0649 +0650 function model = readFieldValue(model, fieldName, value, pos) +0651 if numel(model.(fieldName))<pos-1 +0652 model.(fieldName)(end+1:pos,1) = {''}; +0653 end +0654 model.(fieldName)(pos,1) = {value}; +0655 end +0656 +0657 function [miriams, miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber) +0658 if isempty(key) +0659 key=miriamKey; +0660 else +0661 miriamKey=key; +0662 end +0663 if ~isempty(value) +0664 miriams(entryNumber,1:3) = {pos, key, strip(value)}; +0665 else +0666 entryNumber = entryNumber - 1; +0667 end +0668 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/io/writeYAMLmodel.html b/doc/io/writeYAMLmodel.html index 531ed4c3..39b01959 100644 --- a/doc/io/writeYAMLmodel.html +++ b/doc/io/writeYAMLmodel.html @@ -223,7 +223,7 @@

    SOURCE CODE ^if strcmp(fieldName,'metMiriams') 0166 if ~isempty(model.metMiriams{pos}) -0167 fprintf(fid,[' ' name ': !!omap\n']); +0167 fprintf(fid,' %s: !!omap\n',name); 0168 for i=1:size(model.newMetMiriams,2) 0169 %'i' represents the different miriam names, e.g. 0170 %kegg.compound or chebi @@ -238,7 +238,7 @@

    SOURCE CODE ^elseif strcmp(fieldName,'rxnMiriams') 0181 if ~isempty(model.rxnMiriams{pos}) -0182 fprintf(fid,[' ' name ': !!omap\n']); +0182 fprintf(fid,' %s: !!omap\n',name); 0183 for i=1:size(model.newRxnMiriams,2) 0184 if ~isempty(model.newRxnMiriams{pos,i}) 0185 writeField(model, fid, 'newRxnMiriams', 'txt', pos, [' - ' model.newRxnMiriamNames{i} '_' sprintf('%d',i)], preserveQuotes) @@ -248,7 +248,7 @@

    SOURCE CODE ^elseif strcmp(fieldName,'geneMiriams') 0191 if ~isempty(model.geneMiriams{pos}) -0192 fprintf(fid,[' ' name ': !!omap\n']); +0192 fprintf(fid,' %s: !!omap\n',name); 0193 for i=1:size(model.newGeneMiriams,2) 0194 if ~isempty(model.newGeneMiriams{pos,i}) 0195 writeField(model, fid, 'newGeneMiriams', 'txt', pos, [' - ' model.newGeneMiriamNames{i} '_' sprintf('%d',i)], preserveQuotes) @@ -258,7 +258,7 @@

    SOURCE CODE ^elseif strcmp(fieldName,'compMiriams') 0201 if ~isempty(model.compMiriams{pos}) -0202 fprintf(fid,[' ' name ': !!omap\n']); +0202 fprintf(fid,' %s: !!omap\n',name); 0203 for i=1:size(model.newCompMiriams,2) 0204 if ~isempty(model.newCompMiriams{pos,i}) 0205 writeField(model, fid, 'newCompMiriams', 'txt', pos, [' - ' model.newCompMiriamNames{i} '_' sprintf('%d',i)], preserveQuotes) @@ -268,7 +268,7 @@

    SOURCE CODE ^elseif strcmp(fieldName,'S') 0211 %S: create header & write each metabolite in a new line -0212 fprintf(fid,[' ' name ': !!omap\n']); +0212 fprintf(fid,' %s: !!omap\n',name); 0213 if sum(field(:,pos) ~= 0) > 0 0214 model.mets = model.mets(field(:,pos) ~= 0); 0215 model.coeffs = field(field(:,pos) ~= 0,pos); @@ -282,7 +282,7 @@

    SOURCE CODE ^elseif strcmp(fieldName,'rxnEnzMat') 0225 %S: create header & write each enzyme in a new line -0226 fprintf(fid,[' ' name ': !!omap\n']); +0226 fprintf(fid,' %s: !!omap\n',name); 0227 if sum(field(pos,:) ~= 0) > 0 0228 model.enzymes = model.enzymes(field(pos,:) ~= 0); 0229 model.coeffs = field(pos,field(pos,:) ~= 0); @@ -326,16 +326,16 @@

    SOURCE CODE ^if preserveQuotes 0268 list = ['"' list{1} '"']; 0269 end -0270 fprintf(fid,[' ' name ': ' list '\n']); +0270 fprintf(fid,' %s: %s\n',name,list); 0271 elseif length(list) > 1 || strcmp(fieldName,'subSystems') 0272 if preserveQuotes 0273 for j=1:numel(list) 0274 list{j} = ['"' list{j} '"']; 0275 end 0276 end -0277 fprintf(fid,[' ' name ':\n']); +0277 fprintf(fid,' %s:\n',name); 0278 for i = 1:length(list) -0279 fprintf(fid,[regexprep(name,'(^\s*).*','$1') ' - ' list{i} '\n']); +0279 fprintf(fid,'%s - %s\n',regexprep(name,'(^\s*).*','$1'),list{i}); 0280 end 0281 end 0282 @@ -354,68 +354,65 @@

    SOURCE CODE ^end 0296 end 0297 if ~isempty(value) -0298 fprintf(fid,[' ' name ': ' value '\n']); +0298 fprintf(fid,' %s: %s\n',name,value); 0299 end 0300 end 0301 end -0302 +0302 end 0303 -0304 end -0305 -0306 function writeMetadata(model,fid) -0307 % Writes model metadata to the yaml file. This information will eventually -0308 % be extracted entirely from the model, but for now, many of the entries -0309 % are hard-coded defaults for HumanGEM. -0310 -0311 fprintf(fid, '- metaData:\n'); -0312 fprintf(fid, [' id: "', model.id, '"\n']); -0313 fprintf(fid, [' name: "', model.name, '"\n']); -0314 if isfield(model,'version') -0315 fprintf(fid, [' version: "', model.version, '"\n']); -0316 end -0317 fprintf(fid, [' date: "', datestr(now,29), '"\n']); % 29=YYYY-MM-DD -0318 if isfield(model,'annotation') -0319 if isfield(model.annotation,'defaultLB') -0320 fprintf(fid, [' defaultLB: "', num2str(model.annotation.defaultLB), '"\n']); -0321 end -0322 if isfield(model.annotation,'defaultUB') -0323 fprintf(fid, [' defaultUB: "', num2str(model.annotation.defaultUB), '"\n']); -0324 end -0325 if isfield(model.annotation,'givenName') -0326 fprintf(fid, [' givenName: "', model.annotation.givenName, '"\n']); -0327 end -0328 if isfield(model.annotation,'familyName') -0329 fprintf(fid, [' familyName: "', model.annotation.familyName, '"\n']); -0330 end -0331 if isfield(model.annotation,'authors') -0332 fprintf(fid, [' authors: "', model.annotation.authors, '"\n']); -0333 end -0334 if isfield(model.annotation,'email') -0335 fprintf(fid, [' email: "', model.annotation.email, '"\n']); -0336 end -0337 if isfield(model.annotation,'organization') -0338 fprintf(fid, [' organization: "', model.annotation.organization, '"\n']); -0339 end -0340 if isfield(model.annotation,'taxonomy') -0341 fprintf(fid, [' taxonomy: "', model.annotation.taxonomy, '"\n']); -0342 end -0343 if isfield(model.annotation,'note') -0344 fprintf(fid, [' note: "', model.annotation.note, '"\n']); -0345 end -0346 if isfield(model.annotation,'sourceUrl') -0347 fprintf(fid, [' sourceUrl: "', model.annotation.sourceUrl, '"\n']); -0348 end -0349 end -0350 if isfield(model,'ec') -0351 if model.ec.geckoLight -0352 geckoLight = 'true'; -0353 else -0354 geckoLight = 'false'; -0355 end -0356 fprintf(fid,[' geckoLight: "' geckoLight '"\n']); -0357 end -0358 end -0359 +0304 function writeMetadata(model,fid) +0305 % Writes model metadata to the yaml file. This information will eventually +0306 % be extracted entirely from the model, but for now, many of the entries +0307 % are hard-coded defaults for HumanGEM. +0308 +0309 fprintf(fid, '- metaData:\n'); +0310 fprintf(fid, ' id: "%s"\n', model.id); +0311 fprintf(fid, ' name: "%s"\n',model.name); +0312 if isfield(model,'version') +0313 fprintf(fid, ' version: "%s"\n',model.version); +0314 end +0315 fprintf(fid, ' date: "%s"\n',datestr(now,29)); % 29=YYYY-MM-DD +0316 if isfield(model,'annotation') +0317 if isfield(model.annotation,'defaultLB') +0318 fprintf(fid, ' defaultLB: "%g"\n', model.annotation.defaultLB); +0319 end +0320 if isfield(model.annotation,'defaultUB') +0321 fprintf(fid, ' defaultUB: "%g"\n', model.annotation.defaultUB); +0322 end +0323 if isfield(model.annotation,'givenName') +0324 fprintf(fid, ' givenName: "%s"\n', model.annotation.givenName); +0325 end +0326 if isfield(model.annotation,'familyName') +0327 fprintf(fid, ' familyName: "%s"\n', model.annotation.familyName); +0328 end +0329 if isfield(model.annotation,'authors') +0330 fprintf(fid, ' authors: "%s"\n', model.annotation.authors); +0331 end +0332 if isfield(model.annotation,'email') +0333 fprintf(fid, ' email: "%s"\n', model.annotation.email); +0334 end +0335 if isfield(model.annotation,'organization') +0336 fprintf(fid, ' organization: "%s"\n',model.annotation.organization); +0337 end +0338 if isfield(model.annotation,'taxonomy') +0339 fprintf(fid, ' taxonomy: "%s"\n', model.annotation.taxonomy); +0340 end +0341 if isfield(model.annotation,'note') +0342 fprintf(fid, ' note: "%s"\n', model.annotation.note); +0343 end +0344 if isfield(model.annotation,'sourceUrl') +0345 fprintf(fid, ' sourceUrl: "%s"\n', model.annotation.sourceUrl); +0346 end +0347 end +0348 if isfield(model,'ec') +0349 if model.ec.geckoLight +0350 geckoLight = 'true'; +0351 else +0352 geckoLight = 'false'; +0353 end +0354 fprintf(fid,' geckoLight: "%s"\n',geckoLight); +0355 end +0356 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/doc/struct_conversion/ravenCobraWrapper.html b/doc/struct_conversion/ravenCobraWrapper.html index 6d40c642..9e0a938b 100644 --- a/doc/struct_conversion/ravenCobraWrapper.html +++ b/doc/struct_conversion/ravenCobraWrapper.html @@ -141,369 +141,373 @@

    SOURCE CODE ^'rxnNames','subSystems','rxnNotes',... -0074 'metFormulas','comps','compNames','metCharges','genes',... -0075 'rxnConfidenceScores','rxnGeneMat','metNotes','rev'}; -0076 for i=1:length(optFields) -0077 if isfield(model,optFields{i}) -0078 newModel.(optFields{i})=model.(optFields{i}); -0079 end -0080 end -0081 -0082 % Convert unique fields -0083 if isRaven -0084 fprintf('Converting RAVEN structure to COBRA..\n'); -0085 %Convert from RAVEN to COBRA structure -0086 -0087 %Mandatory COBRA fields -0088 newModel.rxns=model.rxns; -0089 if all(~cellfun(@isempty,regexp(model.mets,'\[[^\]]+\]$'))) -0090 newModel.mets=model.mets; -0091 else -0092 %Check if model has compartment info as "met_c" suffix in all metabolites: -0093 BiGGformat = false(size(model.mets)); -0094 for i=1:numel(model.comps) -0095 compPos=model.metComps==i; -0096 BiGGformat(compPos)=~cellfun(@isempty,regexp(model.mets(compPos),['_' model.comps{i} '$'])); -0097 end -0098 if all(BiGGformat) -0099 newModel.mets=model.mets; -0100 for i=1:numel(model.comps) -0101 newModel.mets=regexprep(newModel.mets,['_' model.comps{i} '$'],['[' model.comps{i} ']']); -0102 end -0103 else -0104 newModel.mets=strcat(model.mets,'[',model.comps(model.metComps),']'); -0105 end -0106 end -0107 -0108 %b, csense, osenseStr, genes, rules are also mandatory, but defined -0109 %later to match the order of fields -0110 -0111 %Optional COBRA fields -0112 if isfield(model,'id') -0113 newModel.modelID=model.id; -0114 end -0115 if isfield(model,'name') -0116 newModel.modelName=model.name; -0117 end -0118 if isfield(model,'eccodes') -0119 newModel.rxnECNumbers=model.eccodes; -0120 end -0121 if isfield(model,'rxnMiriams') -0122 [miriams,extractedMiriamNames]=extractMiriam(model.rxnMiriams); -0123 for i = 1:length(rxnCOBRAfields) -0124 j=ismember(extractedMiriamNames,rxnNamespaces{i}); -0125 if any(j) -0126 eval(['newModel.' rxnCOBRAfields{i} ' = miriams(:,j);']) -0127 end -0128 end -0129 end -0130 if isfield(model,'rxnReferences') % Concatenate model.rxnReferences to those extracted from model.rxnMiriams -0131 if isfield(newModel,'rxnReferences') -0132 newModel.rxnReferences = strcat(newModel.rxnReferences,{'; '},model.rxnReferences); -0133 newModel.rxnReferences = regexprep(newModel.rxnReferences,'^; $',''); -0134 else -0135 newModel.rxnReferences = model.rxnReferences; -0136 end -0137 end -0138 if isfield(model,'metNames') -0139 newModel.metNames=strcat(model.metNames,' [',model.compNames(model.metComps),']'); -0140 end -0141 if isfield(model,'metMiriams') -0142 [miriams,extractedMiriamNames]=extractMiriam(model.metMiriams); -0143 %Shorten miriam names for KEGG and PubChem. These shorter names -0144 %will be used later to concatenate KEGG COMPOUND/GLYCAN and PubChem -0145 %Compound/Substance, into corresponding COBRA model fields -0146 extractedMiriamNames=regexprep(extractedMiriamNames,'^kegg\..+','kegg'); -0147 extractedMiriamNames=regexprep(extractedMiriamNames,'^pubchem\..+','pubchem'); -0148 i=ismember(extractedMiriamNames,'kegg'); -0149 if any(i) % Combine KEGG compounds and glycans -0150 for j=1:length(i) -0151 if i(j) && isfield(newModel,'metKEGGID')~=1 -0152 newModel.metKEGGID=miriams(:,j); -0153 elseif i(j) -0154 newModel.metKEGGID=strcat(newModel.metKEGGID,';',miriams(:,j)); -0155 end -0156 end -0157 newModel.metKEGGID=regexprep(newModel.metKEGGID,'^;|;$',''); -0158 end -0159 i=ismember(extractedMiriamNames,'pubchem'); -0160 if any(i) % Combine Pubchem compounds and substances -0161 for j=1:length(i) -0162 if i(j) && isfield(newModel,'metPubChemID')~=1 -0163 newModel.metPubChemID=miriams(:,j); -0164 elseif i(j) -0165 newModel.metPubChemID=strcat(newModel.metPubChemID,';',miriams(:,j)); -0166 end -0167 end -0168 newModel.metPubChemID=regexprep(newModel.metPubChemID,'^;|;$',''); -0169 end -0170 %All other Miriams can be directly parsed with no modifications: -0171 for i = 1:length(metCOBRAfields) -0172 j=ismember(extractedMiriamNames,metNamespaces{i}); -0173 if any(j) -0174 eval(['newModel.' metCOBRAfields{i} ' = miriams(:,j);']) -0175 end -0176 end -0177 end -0178 if isfield(model,'inchis') -0179 newModel.metInChIString=regexprep(strcat('InChI=', model.inchis),'^InChI=$',''); -0180 end -0181 newModel.b=zeros(numel(model.mets),1); -0182 newModel.csense=repmat('E',size(model.mets)); -0183 if isfield(model,'geneMiriams') -0184 [miriams,extractedMiriamNames]=extractMiriam(model.geneMiriams); -0185 for i = 1:length(geneCOBRAfields) -0186 j=ismember(extractedMiriamNames,geneNamespaces{i}); -0187 if any(j) -0188 eval(['newModel.' geneCOBRAfields{i} ' = miriams(:,j);']) -0189 end -0190 end -0191 end -0192 if isfield(model,'geneShortNames') -0193 newModel.geneNames=model.geneShortNames; -0194 end -0195 if isfield(model,'genes') -0196 newModel.rules=grrulesToRules(model); -0197 else -0198 fprintf('WARNING: no genes detected. The model therefore may not be exportable to SBML file with writeCbModel\n'); -0199 end -0200 newModel.osenseStr='max'; -0201 else -0202 fprintf('Converting COBRA structure to RAVEN..\n'); -0203 %Convert from COBRA to RAVEN structure -0204 -0205 %Mandatory RAVEN fields -0206 newModel.mets=model.mets; -0207 if ~isfield(model,'comps') -0208 model.comps = unique(regexprep(model.mets,'.*\[([^\]]+)\]$','$1')); -0209 end -0210 for i=1:numel(model.comps) -0211 newModel.mets=regexprep(newModel.mets,['\[', model.comps{i}, '\]$'],''); -0212 newModel.mets=regexprep(newModel.mets,['\[', model.compNames{i}, '\]$'],''); +0071 if isfield(model,'c') +0072 newModel.c=model.c; +0073 else +0074 newModel.c=zeros(numel(model.rxns),1); +0075 end +0076 newModel.rxns=model.rxns; +0077 optFields = {'rxnNames','subSystems','rxnNotes',... +0078 'metFormulas','comps','compNames','metCharges','genes',... +0079 'rxnConfidenceScores','rxnGeneMat','metNotes','rev'}; +0080 for i=1:length(optFields) +0081 if isfield(model,optFields{i}) +0082 newModel.(optFields{i})=model.(optFields{i}); +0083 end +0084 end +0085 +0086 % Convert unique fields +0087 if isRaven +0088 fprintf('Converting RAVEN structure to COBRA..\n'); +0089 %Convert from RAVEN to COBRA structure +0090 +0091 %Mandatory COBRA fields +0092 newModel.rxns=model.rxns; +0093 if all(~cellfun(@isempty,regexp(model.mets,'\[[^\]]+\]$'))) +0094 newModel.mets=model.mets; +0095 else +0096 %Check if model has compartment info as "met_c" suffix in all metabolites: +0097 BiGGformat = false(size(model.mets)); +0098 for i=1:numel(model.comps) +0099 compPos=model.metComps==i; +0100 BiGGformat(compPos)=~cellfun(@isempty,regexp(model.mets(compPos),['_' model.comps{i} '$'])); +0101 end +0102 if all(BiGGformat) +0103 newModel.mets=model.mets; +0104 for i=1:numel(model.comps) +0105 newModel.mets=regexprep(newModel.mets,['_' model.comps{i} '$'],['[' model.comps{i} ']']); +0106 end +0107 else +0108 newModel.mets=strcat(model.mets,'[',model.comps(model.metComps),']'); +0109 end +0110 end +0111 +0112 %b, csense, osenseStr, genes, rules are also mandatory, but defined +0113 %later to match the order of fields +0114 +0115 %Optional COBRA fields +0116 if isfield(model,'id') +0117 newModel.modelID=model.id; +0118 end +0119 if isfield(model,'name') +0120 newModel.modelName=model.name; +0121 end +0122 if isfield(model,'eccodes') +0123 newModel.rxnECNumbers=model.eccodes; +0124 end +0125 if isfield(model,'rxnMiriams') +0126 [miriams,extractedMiriamNames]=extractMiriam(model.rxnMiriams); +0127 for i = 1:length(rxnCOBRAfields) +0128 j=ismember(extractedMiriamNames,rxnNamespaces{i}); +0129 if any(j) +0130 eval(['newModel.' rxnCOBRAfields{i} ' = miriams(:,j);']) +0131 end +0132 end +0133 end +0134 if isfield(model,'rxnReferences') % Concatenate model.rxnReferences to those extracted from model.rxnMiriams +0135 if isfield(newModel,'rxnReferences') +0136 newModel.rxnReferences = strcat(newModel.rxnReferences,{'; '},model.rxnReferences); +0137 newModel.rxnReferences = regexprep(newModel.rxnReferences,'^; $',''); +0138 else +0139 newModel.rxnReferences = model.rxnReferences; +0140 end +0141 end +0142 if isfield(model,'metNames') +0143 newModel.metNames=strcat(model.metNames,' [',model.compNames(model.metComps),']'); +0144 end +0145 if isfield(model,'metMiriams') +0146 [miriams,extractedMiriamNames]=extractMiriam(model.metMiriams); +0147 %Shorten miriam names for KEGG and PubChem. These shorter names +0148 %will be used later to concatenate KEGG COMPOUND/GLYCAN and PubChem +0149 %Compound/Substance, into corresponding COBRA model fields +0150 extractedMiriamNames=regexprep(extractedMiriamNames,'^kegg\..+','kegg'); +0151 extractedMiriamNames=regexprep(extractedMiriamNames,'^pubchem\..+','pubchem'); +0152 i=ismember(extractedMiriamNames,'kegg'); +0153 if any(i) % Combine KEGG compounds and glycans +0154 for j=1:length(i) +0155 if i(j) && isfield(newModel,'metKEGGID')~=1 +0156 newModel.metKEGGID=miriams(:,j); +0157 elseif i(j) +0158 newModel.metKEGGID=strcat(newModel.metKEGGID,';',miriams(:,j)); +0159 end +0160 end +0161 newModel.metKEGGID=regexprep(newModel.metKEGGID,'^;|;$',''); +0162 end +0163 i=ismember(extractedMiriamNames,'pubchem'); +0164 if any(i) % Combine Pubchem compounds and substances +0165 for j=1:length(i) +0166 if i(j) && isfield(newModel,'metPubChemID')~=1 +0167 newModel.metPubChemID=miriams(:,j); +0168 elseif i(j) +0169 newModel.metPubChemID=strcat(newModel.metPubChemID,';',miriams(:,j)); +0170 end +0171 end +0172 newModel.metPubChemID=regexprep(newModel.metPubChemID,'^;|;$',''); +0173 end +0174 %All other Miriams can be directly parsed with no modifications: +0175 for i = 1:length(metCOBRAfields) +0176 j=ismember(extractedMiriamNames,metNamespaces{i}); +0177 if any(j) +0178 eval(['newModel.' metCOBRAfields{i} ' = miriams(:,j);']) +0179 end +0180 end +0181 end +0182 if isfield(model,'inchis') +0183 newModel.metInChIString=regexprep(strcat('InChI=', model.inchis),'^InChI=$',''); +0184 end +0185 newModel.b=zeros(numel(model.mets),1); +0186 newModel.csense=repmat('E',size(model.mets)); +0187 if isfield(model,'geneMiriams') +0188 [miriams,extractedMiriamNames]=extractMiriam(model.geneMiriams); +0189 for i = 1:length(geneCOBRAfields) +0190 j=ismember(extractedMiriamNames,geneNamespaces{i}); +0191 if any(j) +0192 eval(['newModel.' geneCOBRAfields{i} ' = miriams(:,j);']) +0193 end +0194 end +0195 end +0196 if isfield(model,'geneShortNames') +0197 newModel.geneNames=model.geneShortNames; +0198 end +0199 if isfield(model,'genes') +0200 newModel.rules=grrulesToRules(model); +0201 else +0202 fprintf('WARNING: no genes detected. The model therefore may not be exportable to SBML file with writeCbModel\n'); +0203 end +0204 newModel.osenseStr='max'; +0205 else +0206 fprintf('Converting COBRA structure to RAVEN..\n'); +0207 %Convert from COBRA to RAVEN structure +0208 +0209 %Mandatory RAVEN fields +0210 newModel.mets=model.mets; +0211 if ~isfield(model,'comps') +0212 model.comps = unique(regexprep(model.mets,'.*\[([^\]]+)\]$','$1')); 0213 end -0214 -0215 %In some cases (e.g. any model that uses BiGG ids as main ids), there -0216 %may be overlapping mets due to removal of compartment info. To avoid -0217 %this, we change compartments from e.g. [c] into _c -0218 if numel(unique(newModel.mets))~=numel(model.mets) -0219 newModel.mets=model.mets; -0220 for i=1:numel(model.comps) -0221 newModel.mets=regexprep(newModel.mets,['\[' model.comps{i} '\]$'],['_' model.comps{i}]); -0222 end -0223 end -0224 %Since COBRA no longer contains rev field it is assumed that rxn is -0225 %reversible if its lower bound is set below zero -0226 if ~isfield(model,'rev') -0227 for i=1:numel(model.rxns) -0228 if model.lb(i)<0 -0229 newModel.rev(i,1)=1; -0230 else -0231 newModel.rev(i,1)=0; -0232 end -0233 end -0234 end -0235 newModel.b=zeros(numel(model.mets),1); -0236 if ~isfield(model,'comps') -0237 %Since 'comps' field is not mandatory in COBRA, it may be required -0238 %to obtain the non-redundant list of comps from metabolite ids, if -0239 %'comps' field is not available -0240 newModel.comps=regexprep(model.mets,'^.+\[',''); -0241 newModel.comps=regexprep(newModel.comps,'\]$',''); -0242 newModel.comps=unique(newModel.comps); -0243 end -0244 -0245 %metComps is also mandatory, but defined later to match the order of -0246 %fields -0247 -0248 %Fields 'name' and 'id' are also considered as mandatory, but -0249 %these are added to the model during exportModel/exportToExcelFormat -0250 %anyway, so there is no point to add this information here +0214 for i=1:numel(model.comps) +0215 newModel.mets=regexprep(newModel.mets,['\[', model.comps{i}, '\]$'],''); +0216 newModel.mets=regexprep(newModel.mets,['\[', model.compNames{i}, '\]$'],''); +0217 end +0218 +0219 %In some cases (e.g. any model that uses BiGG ids as main ids), there +0220 %may be overlapping mets due to removal of compartment info. To avoid +0221 %this, we change compartments from e.g. [c] into _c +0222 if numel(unique(newModel.mets))~=numel(model.mets) +0223 newModel.mets=model.mets; +0224 for i=1:numel(model.comps) +0225 newModel.mets=regexprep(newModel.mets,['\[' model.comps{i} '\]$'],['_' model.comps{i}]); +0226 end +0227 end +0228 %Since COBRA no longer contains rev field it is assumed that rxn is +0229 %reversible if its lower bound is set below zero +0230 if ~isfield(model,'rev') +0231 for i=1:numel(model.rxns) +0232 if model.lb(i)<0 +0233 newModel.rev(i,1)=1; +0234 else +0235 newModel.rev(i,1)=0; +0236 end +0237 end +0238 end +0239 newModel.b=zeros(numel(model.mets),1); +0240 if ~isfield(model,'comps') +0241 %Since 'comps' field is not mandatory in COBRA, it may be required +0242 %to obtain the non-redundant list of comps from metabolite ids, if +0243 %'comps' field is not available +0244 newModel.comps=regexprep(model.mets,'^.+\[',''); +0245 newModel.comps=regexprep(newModel.comps,'\]$',''); +0246 newModel.comps=unique(newModel.comps); +0247 end +0248 +0249 %metComps is also mandatory, but defined later to match the order of +0250 %fields 0251 -0252 %Optional RAVEN fields -0253 if isfield(model,'modelID') -0254 newModel.id=model.modelID; -0255 end -0256 if isfield(model,'modelName') -0257 newModel.name=model.modelName; -0258 end -0259 if isfield(model,'rules') -0260 model.grRules = rulesTogrrules(model); -0261 [grRules,rxnGeneMat] = standardizeGrRules(model,true); -0262 newModel.grRules = grRules; -0263 newModel.rxnGeneMat = rxnGeneMat; -0264 end -0265 if isfield(model,'rxnECNumbers') -0266 newModel.eccodes=regexprep(model.rxnECNumbers,'EC|EC:',''); -0267 end -0268 if any(isfield(model,rxnCOBRAfields)) -0269 for i=1:numel(model.rxns) -0270 counter=1; -0271 newModel.rxnMiriams{i,1}=[]; -0272 if isfield(model,'rxnReferences') -0273 if ~isempty(model.rxnReferences{i}) -0274 pmids = model.rxnReferences{i}; -0275 pmids = strsplit(pmids,'; '); -0276 nonPmids = cellfun(@isempty,regexp(pmids,'^\d+$','match','once')); -0277 if any(nonPmids) %Not a pubmed id, keep in rxnReferences instead -0278 newModel.rxnReferences{i,1} = strjoin(pmids(nonPmids),', '); -0279 pmids(nonPmids)=[]; -0280 end -0281 for j = 1:length(pmids) -0282 newModel.rxnMiriams{i,1}.name{counter,1} = 'pubmed'; -0283 newModel.rxnMiriams{i,1}.value{counter,1} = pmids{j}; -0284 counter=counter+1; -0285 end -0286 end -0287 end -0288 for j = 2:length(rxnCOBRAfields) %Start from 2, as 1 is rxnReferences -0289 if isfield(model,rxnCOBRAfields{j}) -0290 rxnAnnotation = eval(['model.' rxnCOBRAfields{j} '{i}']); -0291 if ~isempty(rxnAnnotation) -0292 rxnAnnotation = strtrim(strsplit(rxnAnnotation,';')); -0293 for a=1:length(rxnAnnotation) -0294 newModel.rxnMiriams{i,1}.name{counter,1} = rxnNamespaces{j}; -0295 newModel.rxnMiriams{i,1}.value{counter,1} = rxnAnnotation{a}; -0296 counter=counter+1; -0297 end -0298 end -0299 end -0300 end -0301 end -0302 end -0303 if isfield(newModel,'rxnReferences') -0304 emptyEntry = cellfun(@isempty,newModel.rxnReferences); -0305 newModel.rxnReferences(emptyEntry)={''}; +0252 %Fields 'name' and 'id' are also considered as mandatory, but +0253 %these are added to the model during exportModel/exportToExcelFormat +0254 %anyway, so there is no point to add this information here +0255 +0256 %Optional RAVEN fields +0257 if isfield(model,'modelID') +0258 newModel.id=model.modelID; +0259 end +0260 if isfield(model,'modelName') +0261 newModel.name=model.modelName; +0262 end +0263 if isfield(model,'rules') +0264 model.grRules = rulesTogrrules(model); +0265 [grRules,rxnGeneMat] = standardizeGrRules(model,true); +0266 newModel.grRules = grRules; +0267 newModel.rxnGeneMat = rxnGeneMat; +0268 end +0269 if isfield(model,'rxnECNumbers') +0270 newModel.eccodes=regexprep(model.rxnECNumbers,'EC|EC:',''); +0271 end +0272 if any(isfield(model,rxnCOBRAfields)) +0273 for i=1:numel(model.rxns) +0274 counter=1; +0275 newModel.rxnMiriams{i,1}=[]; +0276 if isfield(model,'rxnReferences') +0277 if ~isempty(model.rxnReferences{i}) +0278 pmids = model.rxnReferences{i}; +0279 pmids = strsplit(pmids,'; '); +0280 nonPmids = cellfun(@isempty,regexp(pmids,'^\d+$','match','once')); +0281 if any(nonPmids) %Not a pubmed id, keep in rxnReferences instead +0282 newModel.rxnReferences{i,1} = strjoin(pmids(nonPmids),', '); +0283 pmids(nonPmids)=[]; +0284 end +0285 for j = 1:length(pmids) +0286 newModel.rxnMiriams{i,1}.name{counter,1} = 'pubmed'; +0287 newModel.rxnMiriams{i,1}.value{counter,1} = pmids{j}; +0288 counter=counter+1; +0289 end +0290 end +0291 end +0292 for j = 2:length(rxnCOBRAfields) %Start from 2, as 1 is rxnReferences +0293 if isfield(model,rxnCOBRAfields{j}) +0294 rxnAnnotation = eval(['model.' rxnCOBRAfields{j} '{i}']); +0295 if ~isempty(rxnAnnotation) +0296 rxnAnnotation = strtrim(strsplit(rxnAnnotation,';')); +0297 for a=1:length(rxnAnnotation) +0298 newModel.rxnMiriams{i,1}.name{counter,1} = rxnNamespaces{j}; +0299 newModel.rxnMiriams{i,1}.value{counter,1} = rxnAnnotation{a}; +0300 counter=counter+1; +0301 end +0302 end +0303 end +0304 end +0305 end 0306 end -0307 if any(isfield(model,geneCOBRAfields)) -0308 for i=1:numel(model.genes) -0309 counter=1; -0310 newModel.geneMiriams{i,1}=[]; -0311 for j = 1:length(geneCOBRAfields) -0312 if isfield(model,geneCOBRAfields{j}) -0313 geneAnnotation = eval(['model.' geneCOBRAfields{j} '{i}']); -0314 if ~isempty(geneAnnotation) -0315 geneAnnotation = strtrim(strsplit(geneAnnotation,';')); -0316 for a=1:length(geneAnnotation) -0317 newModel.geneMiriams{i,1}.name{counter,1} = geneNamespaces{j}; -0318 newModel.geneMiriams{i,1}.value{counter,1} = geneAnnotation{a}; -0319 counter=counter+1; -0320 end -0321 end -0322 end -0323 end -0324 end -0325 end -0326 if isfield(model,'geneNames') -0327 newModel.geneShortNames=model.geneNames; -0328 end -0329 newModel.metNames=model.metNames; -0330 for i=1:numel(model.comps) -0331 newModel.metNames=regexprep(newModel.metNames,['\[', model.comps{i}, '\]$'],''); -0332 newModel.metNames=regexprep(newModel.metNames,['\[', model.compNames{i}, '\]$'],''); -0333 end -0334 newModel.metNames=deblank(newModel.metNames); -0335 newModel.metComps=regexprep(model.mets,'^.+\[',''); -0336 newModel.metComps=regexprep(newModel.metComps,'\]$',''); -0337 [~, newModel.metComps]=ismember(newModel.metComps,newModel.comps); -0338 if isfield(model,'metInChIString') -0339 newModel.inchis=regexprep(model.metInChIString,'^InChI=',''); -0340 end -0341 printWarning=false; -0342 if any(isfield(model,[metCOBRAfields;'metKEGGID';'metPubChemID'])) -0343 for i=1:numel(model.mets) -0344 counter=1; -0345 newModel.metMiriams{i,1}=[]; -0346 if isfield(model,'metKEGGID') -0347 if ~isempty(model.metKEGGID{i}) -0348 if strcmp(model.metKEGGID{i}(1),'C') -0349 newModel.metMiriams{i,1}.name{counter,1} = 'kegg.compound'; -0350 newModel.metMiriams{i,1}.value{counter,1} = model.metKEGGID{i}; -0351 counter=counter+1; -0352 elseif strcmp(model.metKEGGID{i}(1),'G') -0353 newModel.metMiriams{i,1}.name{counter,1} = 'kegg.glycan'; +0307 if isfield(newModel,'rxnReferences') +0308 emptyEntry = cellfun(@isempty,newModel.rxnReferences); +0309 newModel.rxnReferences(emptyEntry)={''}; +0310 end +0311 if any(isfield(model,geneCOBRAfields)) +0312 for i=1:numel(model.genes) +0313 counter=1; +0314 newModel.geneMiriams{i,1}=[]; +0315 for j = 1:length(geneCOBRAfields) +0316 if isfield(model,geneCOBRAfields{j}) +0317 geneAnnotation = eval(['model.' geneCOBRAfields{j} '{i}']); +0318 if ~isempty(geneAnnotation) +0319 geneAnnotation = strtrim(strsplit(geneAnnotation,';')); +0320 for a=1:length(geneAnnotation) +0321 newModel.geneMiriams{i,1}.name{counter,1} = geneNamespaces{j}; +0322 newModel.geneMiriams{i,1}.value{counter,1} = geneAnnotation{a}; +0323 counter=counter+1; +0324 end +0325 end +0326 end +0327 end +0328 end +0329 end +0330 if isfield(model,'geneNames') +0331 newModel.geneShortNames=model.geneNames; +0332 end +0333 newModel.metNames=model.metNames; +0334 for i=1:numel(model.comps) +0335 newModel.metNames=regexprep(newModel.metNames,['\[', model.comps{i}, '\]$'],''); +0336 newModel.metNames=regexprep(newModel.metNames,['\[', model.compNames{i}, '\]$'],''); +0337 end +0338 newModel.metNames=deblank(newModel.metNames); +0339 newModel.metComps=regexprep(model.mets,'^.+\[',''); +0340 newModel.metComps=regexprep(newModel.metComps,'\]$',''); +0341 [~, newModel.metComps]=ismember(newModel.metComps,newModel.comps); +0342 if isfield(model,'metInChIString') +0343 newModel.inchis=regexprep(model.metInChIString,'^InChI=',''); +0344 end +0345 printWarning=false; +0346 if any(isfield(model,[metCOBRAfields;'metKEGGID';'metPubChemID'])) +0347 for i=1:numel(model.mets) +0348 counter=1; +0349 newModel.metMiriams{i,1}=[]; +0350 if isfield(model,'metKEGGID') +0351 if ~isempty(model.metKEGGID{i}) +0352 if strcmp(model.metKEGGID{i}(1),'C') +0353 newModel.metMiriams{i,1}.name{counter,1} = 'kegg.compound'; 0354 newModel.metMiriams{i,1}.value{counter,1} = model.metKEGGID{i}; 0355 counter=counter+1; -0356 end -0357 end -0358 end -0359 if isfield(model,'metPubChemID') -0360 if ~isempty(model.metPubChemID{i}) -0361 if length(model.metPubChemID{i})>3 && strcmp(model.metPubChemID{i}(1:4),'CID:') -0362 newModel.metMiriams{i,1}.name{counter,1} = 'pubchem.compound'; -0363 newModel.metMiriams{i,1}.value{counter,1} = model.metPubChemID{i}; -0364 counter=counter+1; -0365 elseif length(model.metPubChemID{i})>3 && strcmp(model.metPubChemID{i}(1:4),'SID:') -0366 newModel.metMiriams{i,1}.name{counter,1} = 'pubchem.substance'; +0356 elseif strcmp(model.metKEGGID{i}(1),'G') +0357 newModel.metMiriams{i,1}.name{counter,1} = 'kegg.glycan'; +0358 newModel.metMiriams{i,1}.value{counter,1} = model.metKEGGID{i}; +0359 counter=counter+1; +0360 end +0361 end +0362 end +0363 if isfield(model,'metPubChemID') +0364 if ~isempty(model.metPubChemID{i}) +0365 if length(model.metPubChemID{i})>3 && strcmp(model.metPubChemID{i}(1:4),'CID:') +0366 newModel.metMiriams{i,1}.name{counter,1} = 'pubchem.compound'; 0367 newModel.metMiriams{i,1}.value{counter,1} = model.metPubChemID{i}; 0368 counter=counter+1; -0369 else -0370 newModel.metMiriams{i,1}.name{counter,1} = 'pubchem.compound'; +0369 elseif length(model.metPubChemID{i})>3 && strcmp(model.metPubChemID{i}(1:4),'SID:') +0370 newModel.metMiriams{i,1}.name{counter,1} = 'pubchem.substance'; 0371 newModel.metMiriams{i,1}.value{counter,1} = model.metPubChemID{i}; 0372 counter=counter+1; -0373 printWarning=true; -0374 end -0375 end -0376 end -0377 for j = 1:length(metCOBRAfields) -0378 if isfield(model,metCOBRAfields{j}) -0379 metAnnotation = eval(['model.' metCOBRAfields{j} '{i}']); -0380 if ~isempty(metAnnotation) -0381 metAnnotation = strtrim(strsplit(metAnnotation,';')); -0382 for a=1:length(metAnnotation) -0383 newModel.metMiriams{i,1}.name{counter,1} = metNamespaces{j}; -0384 newModel.metMiriams{i,1}.value{counter,1} = metAnnotation{a}; -0385 counter=counter+1; -0386 end -0387 end -0388 end -0389 end -0390 end -0391 end -0392 if printWarning -0393 fprintf('Could not determine whether PubChemIDs are compounds (CID)\n or substances (SID). All annotated PubChemIDs will therefore \n be assigned as compounds (CID).\n'); -0394 end -0395 end -0396 -0397 % Order fields -0398 modelNew=standardizeModelFieldOrder(newModel); % Corrects for both RAVEN and COBRA models +0373 else +0374 newModel.metMiriams{i,1}.name{counter,1} = 'pubchem.compound'; +0375 newModel.metMiriams{i,1}.value{counter,1} = model.metPubChemID{i}; +0376 counter=counter+1; +0377 printWarning=true; +0378 end +0379 end +0380 end +0381 for j = 1:length(metCOBRAfields) +0382 if isfield(model,metCOBRAfields{j}) +0383 metAnnotation = eval(['model.' metCOBRAfields{j} '{i}']); +0384 if ~isempty(metAnnotation) +0385 metAnnotation = strtrim(strsplit(metAnnotation,';')); +0386 for a=1:length(metAnnotation) +0387 newModel.metMiriams{i,1}.name{counter,1} = metNamespaces{j}; +0388 newModel.metMiriams{i,1}.value{counter,1} = metAnnotation{a}; +0389 counter=counter+1; +0390 end +0391 end +0392 end +0393 end +0394 end +0395 end +0396 if printWarning +0397 fprintf('Could not determine whether PubChemIDs are compounds (CID)\n or substances (SID). All annotated PubChemIDs will therefore \n be assigned as compounds (CID).\n'); +0398 end 0399 end 0400 -0401 function rules=grrulesToRules(model) -0402 %This function just takes grRules, changes all gene names to -0403 %'x(geneNumber)' and also changes 'or' and 'and' relations to corresponding -0404 %symbols -0405 replacingGenes=cell([size(model.genes,1) 1]); -0406 for i=1:numel(replacingGenes) -0407 replacingGenes{i}=strcat('x(',num2str(i),')'); -0408 end -0409 rules = strcat({' '},model.grRules,{' '}); -0410 for i=1:length(model.genes) -0411 rules=regexprep(rules,[' ' model.genes{i} ' '],[' ' replacingGenes{i} ' ']); -0412 rules=regexprep(rules,['(' model.genes{i} ' '],['(' replacingGenes{i} ' ']); -0413 rules=regexprep(rules,[' ' model.genes{i} ')'],[' ' replacingGenes{i} ')']); -0414 end -0415 rules=regexprep(rules,' and ',' & '); -0416 rules=regexprep(rules,' or ',' | '); -0417 rules=strtrim(rules); +0401 % Order fields +0402 modelNew=standardizeModelFieldOrder(newModel); % Corrects for both RAVEN and COBRA models +0403 end +0404 +0405 function rules=grrulesToRules(model) +0406 %This function just takes grRules, changes all gene names to +0407 %'x(geneNumber)' and also changes 'or' and 'and' relations to corresponding +0408 %symbols +0409 replacingGenes=cell([size(model.genes,1) 1]); +0410 for i=1:numel(replacingGenes) +0411 replacingGenes{i}=strcat('x(',num2str(i),')'); +0412 end +0413 rules = strcat({' '},model.grRules,{' '}); +0414 for i=1:length(model.genes) +0415 rules=regexprep(rules,[' ' model.genes{i} ' '],[' ' replacingGenes{i} ' ']); +0416 rules=regexprep(rules,['(' model.genes{i} ' '],['(' replacingGenes{i} ' ']); +0417 rules=regexprep(rules,[' ' model.genes{i} ')'],[' ' replacingGenes{i} ')']); 0418 end -0419 -0420 function grRules=rulesTogrrules(model) -0421 %This function takes rules, replaces &/| for and/or, replaces the x(i) -0422 %format with the actual gene ID, and takes out extra whitespace and -0423 %redundant parenthesis introduced by COBRA, to create grRules. -0424 grRules = strrep(model.rules,'&','and'); -0425 grRules = strrep(grRules,'|','or'); -0426 for i = 1:length(model.genes) -0427 grRules = strrep(grRules,['x(' num2str(i) ')'],model.genes{i}); -0428 end -0429 grRules = strrep(grRules,'( ','('); -0430 grRules = strrep(grRules,' )',')'); -0431 grRules = regexprep(grRules,'^(',''); %rules that start with a "(" -0432 grRules = regexprep(grRules,')$',''); %rules that end with a ")" -0433 end +0419 rules=regexprep(rules,' and ',' & '); +0420 rules=regexprep(rules,' or ',' | '); +0421 rules=strtrim(rules); +0422 end +0423 +0424 function grRules=rulesTogrrules(model) +0425 %This function takes rules, replaces &/| for and/or, replaces the x(i) +0426 %format with the actual gene ID, and takes out extra whitespace and +0427 %redundant parenthesis introduced by COBRA, to create grRules. +0428 grRules = strrep(model.rules,'&','and'); +0429 grRules = strrep(grRules,'|','or'); +0430 for i = 1:length(model.genes) +0431 grRules = strrep(grRules,['x(' num2str(i) ')'],model.genes{i}); +0432 end +0433 grRules = strrep(grRules,'( ','('); +0434 grRules = strrep(grRules,' )',')'); +0435 grRules = regexprep(grRules,'^(',''); %rules that start with a "(" +0436 grRules = regexprep(grRules,')$',''); %rules that end with a ")" +0437 end
    Generated by m2html © 2005
    \ No newline at end of file diff --git a/external/kegg/getPhylDist.m b/external/kegg/getPhylDist.m index f9f5ae0f..02b8f18d 100755 --- a/external/kegg/getPhylDist.m +++ b/external/kegg/getPhylDist.m @@ -48,6 +48,7 @@ fid = fopen(fullfile(keggPath,'taxonomy'), 'r'); phylDistStruct.ids={}; + phylDistStruct.names={}; %Keeps the categories for each organism orgCat={}; @@ -93,6 +94,7 @@ %Should always exist phylDistStruct.ids{orgCounter}=tline(sPos(1)+1:sPos(2)-1); + phylDistStruct.names{orgCounter}=tline(sPos(3)+1:end); orgCat{orgCounter}=currentCat; end end diff --git a/external/kegg/keggPhylDist.mat b/external/kegg/keggPhylDist.mat index 36f1dd5e..9cdde3b9 100644 Binary files a/external/kegg/keggPhylDist.mat and b/external/kegg/keggPhylDist.mat differ diff --git a/external/metacyc/getMetaCycModelForOrganism.m b/external/metacyc/getMetaCycModelForOrganism.m index aa23ded1..ac954864 100755 --- a/external/metacyc/getMetaCycModelForOrganism.m +++ b/external/metacyc/getMetaCycModelForOrganism.m @@ -1,5 +1,5 @@ -function model=getMetaCycModelForOrganism(organismID,fastaFile,... - keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond) +function [model, metacycBlastStruct] = getMetaCycModelForOrganism(organismID,fastaFile,... + keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond,metacycBlastStruct) % getMetaCycModelForOrganism % Reconstructs a genome-scale metabolic model based on protein homology to the % MetaCyc pathway database @@ -23,12 +23,19 @@ % minPositives minimum Positives values of BLASTp search (opt, default 45 %) % useDiamond use DIAMOND alignment tools to perform homology search % if true, otherwise the BLASTP is used (opt, default true) +% metacycBlastStruct provided an earlier generated metacycBlastStruct +% from getMetaCycModelForOrganism (opt, default empty). +% Useful when trying different cutoffs for minScore +% and minPositives without having to regenerate the +% blastStructure each time. % % Output: % model a model structure for the query organism +% metacycBlastStruct result from getBlast or getDiamond, before the +% minScore and minPositives cutoffs are applied % -% Usage: model=getMetaCycModelForOrganism(organismID,fastaFile,... -% keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond) +% Usage: [model, metacycBlastStruct] = getMetaCycModelForOrganism(organismID,fastaFile,... +% keepTransportRxns,keepUnbalanced,keepUndetermined,minScore,minPositives,useDiamond,metacycBlastStruct) organismID=char(organismID); if nargin<2 @@ -55,6 +62,9 @@ if nargin<8 useDiamond=true; end +if nargin<9 + metacycBlastStruct=[]; +end if isempty(fastaFile) error('*** The query FASTA filename cannot be empty! ***'); @@ -88,15 +98,19 @@ ravenPath=findRAVENroot(); %Generate blast strcture by either DIAMOND or BLASTP -if useDiamond - blastStruc=getDiamond(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa')); +if isempty(metacycBlastStruct) + if useDiamond + blastStruc=getDiamond(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa')); + else + blastStruc=getBlast(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa')); + end + %Only look the query + blastStructure=blastStruc(2); + metacycBlastStruct=blastStructure; else - blastStruc=getBlast(organismID,fastaFile,{'MetaCyc'},fullfile(ravenPath,'external','metacyc','protseq.fsa')); + blastStructure=metacycBlastStruct; end -%Only look the query -blastStructure=blastStruc(2); - %Remove all blast hits that are below the cutoffs indexes=blastStructure.bitscore>=minScore & blastStructure.ppos>=minPositives; blastStructure.toGenes(~indexes)=[]; diff --git a/installation/checkFunctionUniqueness.m b/installation/checkFunctionUniqueness.m index 87019777..d4f9da9f 100755 --- a/installation/checkFunctionUniqueness.m +++ b/installation/checkFunctionUniqueness.m @@ -23,7 +23,12 @@ ravenFunctions={temp_res1.name,temp_res2.name}'; else - altDirs = strsplit(altDirs,';')'; + if ispc + splitChar = ';'; + else + splitChar = ':'; + end + altDirs = strsplit(altDirs,splitChar)'; ravenFunctions={}; for i=1:numel(altDirs) temp_res=dir(fullfile(char(altDirs(i)),'*.m')); diff --git a/installation/checkInstallation.m b/installation/checkInstallation.m index aa5223f1..d97dff16 100755 --- a/installation/checkInstallation.m +++ b/installation/checkInstallation.m @@ -339,3 +339,15 @@ end end end + +function printOrange(stringToPrint) +% printOrange +% Duplicate of RAVEN/core/printOrange is also kept here, as this function +% should be able to run before adding RAVEN to the MATLAB path. +try useDesktop = usejava('desktop'); catch, useDesktop = false; end +if useDesktop + fprintf(['[\b' stringToPrint,']\b']) +else + fprintf(stringToPrint) +end +end diff --git a/io/importModel.m b/io/importModel.m index 1c6fe1c8..f0b23785 100755 --- a/io/importModel.m +++ b/io/importModel.m @@ -525,6 +525,12 @@ lb=regexprep(lb,parameter.name(n),num2str(parameter.value{n})); ub=regexprep(ub,parameter.name(n),num2str(parameter.value{n})); end + if isempty(lb) + lb='-Inf'; + end + if isempty(ub) + ub='Inf'; + end reactionLB(counter)=str2num(lb); reactionUB(counter)=str2num(ub); %The order of these parameters should not be hard coded diff --git a/io/readYAMLmodel.m b/io/readYAMLmodel.m index a56420e3..a32df0ff 100644 --- a/io/readYAMLmodel.m +++ b/io/readYAMLmodel.m @@ -578,14 +578,15 @@ % end % end -% Make rxnGeneMat fields -[genes, rxnGeneMat] = getGenesFromGrRules(model.grRules, model.genes); -if isequal(sort(genes), sort(model.genes)) - model.rxnGeneMat = rxnGeneMat; - model.genes = genes; -else - error('The gene list and grRules are inconsistent.'); -end +% Make rxnGeneMat fields and map to the existing model.genes field +[genes, rxnGeneMat] = getGenesFromGrRules(model.grRules); +model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes)); +[~,geneOrder] = ismember(genes,model.genes); +if any(geneOrder == 0) + error(['The grRules includes the following gene(s), that are not in '... + 'the list of model genes: ', genes{~geneOrder}]) +end +model.rxnGeneMat(:,geneOrder) = rxnGeneMat; % Finalize GECKO model if isGECKO diff --git a/io/writeYAMLmodel.m b/io/writeYAMLmodel.m index b1ac8b34..e50eb4eb 100644 --- a/io/writeYAMLmodel.m +++ b/io/writeYAMLmodel.m @@ -164,7 +164,7 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) if strcmp(fieldName,'metMiriams') if ~isempty(model.metMiriams{pos}) - fprintf(fid,[' ' name ': !!omap\n']); + fprintf(fid,' %s: !!omap\n',name); for i=1:size(model.newMetMiriams,2) %'i' represents the different miriam names, e.g. %kegg.compound or chebi @@ -179,7 +179,7 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) elseif strcmp(fieldName,'rxnMiriams') if ~isempty(model.rxnMiriams{pos}) - fprintf(fid,[' ' name ': !!omap\n']); + fprintf(fid,' %s: !!omap\n',name); for i=1:size(model.newRxnMiriams,2) if ~isempty(model.newRxnMiriams{pos,i}) writeField(model, fid, 'newRxnMiriams', 'txt', pos, [' - ' model.newRxnMiriamNames{i} '_' sprintf('%d',i)], preserveQuotes) @@ -189,7 +189,7 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) elseif strcmp(fieldName,'geneMiriams') if ~isempty(model.geneMiriams{pos}) - fprintf(fid,[' ' name ': !!omap\n']); + fprintf(fid,' %s: !!omap\n',name); for i=1:size(model.newGeneMiriams,2) if ~isempty(model.newGeneMiriams{pos,i}) writeField(model, fid, 'newGeneMiriams', 'txt', pos, [' - ' model.newGeneMiriamNames{i} '_' sprintf('%d',i)], preserveQuotes) @@ -199,7 +199,7 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) elseif strcmp(fieldName,'compMiriams') if ~isempty(model.compMiriams{pos}) - fprintf(fid,[' ' name ': !!omap\n']); + fprintf(fid,' %s: !!omap\n',name); for i=1:size(model.newCompMiriams,2) if ~isempty(model.newCompMiriams{pos,i}) writeField(model, fid, 'newCompMiriams', 'txt', pos, [' - ' model.newCompMiriamNames{i} '_' sprintf('%d',i)], preserveQuotes) @@ -209,7 +209,7 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) elseif strcmp(fieldName,'S') %S: create header & write each metabolite in a new line - fprintf(fid,[' ' name ': !!omap\n']); + fprintf(fid,' %s: !!omap\n',name); if sum(field(:,pos) ~= 0) > 0 model.mets = model.mets(field(:,pos) ~= 0); model.coeffs = field(field(:,pos) ~= 0,pos); @@ -223,7 +223,7 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) elseif strcmp(fieldName,'rxnEnzMat') %S: create header & write each enzyme in a new line - fprintf(fid,[' ' name ': !!omap\n']); + fprintf(fid,' %s: !!omap\n',name); if sum(field(pos,:) ~= 0) > 0 model.enzymes = model.enzymes(field(pos,:) ~= 0); model.coeffs = field(pos,field(pos,:) ~= 0); @@ -267,16 +267,16 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) if preserveQuotes list = ['"' list{1} '"']; end - fprintf(fid,[' ' name ': ' list '\n']); + fprintf(fid,' %s: %s\n',name,list); elseif length(list) > 1 || strcmp(fieldName,'subSystems') if preserveQuotes for j=1:numel(list) list{j} = ['"' list{j} '"']; end end - fprintf(fid,[' ' name ':\n']); + fprintf(fid,' %s:\n',name); for i = 1:length(list) - fprintf(fid,[regexprep(name,'(^\s*).*','$1') ' - ' list{i} '\n']); + fprintf(fid,'%s - %s\n',regexprep(name,'(^\s*).*','$1'),list{i}); end end @@ -295,12 +295,10 @@ function writeField(model,fid,fieldName,type,pos,name,preserveQuotes) end end if ~isempty(value) - fprintf(fid,[' ' name ': ' value '\n']); + fprintf(fid,' %s: %s\n',name,value); end end end - - end function writeMetadata(model,fid) @@ -309,42 +307,42 @@ function writeMetadata(model,fid) % are hard-coded defaults for HumanGEM. fprintf(fid, '- metaData:\n'); -fprintf(fid, [' id: "', model.id, '"\n']); -fprintf(fid, [' name: "', model.name, '"\n']); +fprintf(fid, ' id: "%s"\n', model.id); +fprintf(fid, ' name: "%s"\n',model.name); if isfield(model,'version') - fprintf(fid, [' version: "', model.version, '"\n']); + fprintf(fid, ' version: "%s"\n',model.version); end -fprintf(fid, [' date: "', datestr(now,29), '"\n']); % 29=YYYY-MM-DD +fprintf(fid, ' date: "%s"\n',datestr(now,29)); % 29=YYYY-MM-DD if isfield(model,'annotation') if isfield(model.annotation,'defaultLB') - fprintf(fid, [' defaultLB: "', num2str(model.annotation.defaultLB), '"\n']); + fprintf(fid, ' defaultLB: "%g"\n', model.annotation.defaultLB); end if isfield(model.annotation,'defaultUB') - fprintf(fid, [' defaultUB: "', num2str(model.annotation.defaultUB), '"\n']); + fprintf(fid, ' defaultUB: "%g"\n', model.annotation.defaultUB); end if isfield(model.annotation,'givenName') - fprintf(fid, [' givenName: "', model.annotation.givenName, '"\n']); + fprintf(fid, ' givenName: "%s"\n', model.annotation.givenName); end if isfield(model.annotation,'familyName') - fprintf(fid, [' familyName: "', model.annotation.familyName, '"\n']); + fprintf(fid, ' familyName: "%s"\n', model.annotation.familyName); end if isfield(model.annotation,'authors') - fprintf(fid, [' authors: "', model.annotation.authors, '"\n']); + fprintf(fid, ' authors: "%s"\n', model.annotation.authors); end if isfield(model.annotation,'email') - fprintf(fid, [' email: "', model.annotation.email, '"\n']); + fprintf(fid, ' email: "%s"\n', model.annotation.email); end if isfield(model.annotation,'organization') - fprintf(fid, [' organization: "', model.annotation.organization, '"\n']); + fprintf(fid, ' organization: "%s"\n',model.annotation.organization); end if isfield(model.annotation,'taxonomy') - fprintf(fid, [' taxonomy: "', model.annotation.taxonomy, '"\n']); + fprintf(fid, ' taxonomy: "%s"\n', model.annotation.taxonomy); end if isfield(model.annotation,'note') - fprintf(fid, [' note: "', model.annotation.note, '"\n']); + fprintf(fid, ' note: "%s"\n', model.annotation.note); end if isfield(model.annotation,'sourceUrl') - fprintf(fid, [' sourceUrl: "', model.annotation.sourceUrl, '"\n']); + fprintf(fid, ' sourceUrl: "%s"\n', model.annotation.sourceUrl); end end if isfield(model,'ec') @@ -353,7 +351,6 @@ function writeMetadata(model,fid) else geckoLight = 'false'; end - fprintf(fid,[' geckoLight: "' geckoLight '"\n']); + fprintf(fid,' geckoLight: "%s"\n',geckoLight); end end - diff --git a/struct_conversion/ravenCobraWrapper.m b/struct_conversion/ravenCobraWrapper.m index 1d321139..196c6d67 100755 --- a/struct_conversion/ravenCobraWrapper.m +++ b/struct_conversion/ravenCobraWrapper.m @@ -68,7 +68,11 @@ newModel.S=model.S; newModel.lb=model.lb; newModel.ub=model.ub; -newModel.c=model.c; +if isfield(model,'c') + newModel.c=model.c; +else + newModel.c=zeros(numel(model.rxns),1); +end newModel.rxns=model.rxns; optFields = {'rxnNames','subSystems','rxnNotes',... 'metFormulas','comps','compNames','metCharges','genes',...