readYAMLmodel

 readYAMLmodel
   Reads a yaml file matching (roughly) the cobrapy yaml structure

   Input:
   fileName    a model file in yaml file format. A dialog window will open
               if no file name is specified.
   verbose     set as true to monitor progress (optional, default false)

   Output:
   model       a model structure

 Usage: model = readYAMLmodel(fileName, verbose)

0001 function model=readYAMLmodel(fileName, verbose)
0002 % readYAMLmodel
0003 %   Reads a yaml file matching (roughly) the cobrapy yaml structure
0004 %
0005 %   Input:
0006 %   fileName    a model file in yaml file format. A dialog window will open
0007 %               if no file name is specified.
0008 %   verbose     set as true to monitor progress (optional, default false)
0009 %
0010 %   Output:
0011 %   model       a model structure
0012 %
0013 % Usage: model = readYAMLmodel(fileName, verbose)
0014 if nargin<1 || isempty(fileName)
0015     [fileName, pathName] = uigetfile({'*.yml;*.yaml'}, 'Please select the model file');
0016     if fileName == 0
0017         error('You should select a model file')
0018     else
0019         fileName = fullfile(pathName,fileName);
0020     end
0021 end
0022 if nargin < 2
0023     verbose = false;
0024 end
0025 
0026 if ~isfile(fileName)
0027     error('Yaml file %s cannot be found', string(fileName));
0028 end
0029 
0030 if verLessThan('matlab','9.9') %readlines introduced 2020b
0031     fid=fopen(fileName);
0032     line_raw=cell(1000000,1);
0033     i=1;
0034     while ~feof(fid)
0035         line_raw{i}=fgetl(fid);
0036         i=i+1;
0037     end
0038     line_raw(i:end)=[];
0039     line_raw=string(line_raw);
0040     fclose(fid);
0041 else
0042     line_raw=readlines(fileName);
0043 end
0044 % If entry is broken of multiple lines, concatenate. Assumes at least 6
0045 % leading spaces to avoid metaData to be concatenated.
0046 newLine=regexp(line_raw,'^ {6,}([\w\(\)].*)','tokenExtents');
0047 brokenLine=find(~cellfun('isempty',newLine));
0048 for i=flip(1:numel(brokenLine))
0049     extraLine = char(line_raw(brokenLine(i)));
0050     extraLine = extraLine(newLine{brokenLine(i)}{1}(1):end);
0051     line_raw{brokenLine(i)-1} = strjoin({line_raw{brokenLine(i)-1},extraLine},' ');
0052 end
0053 line_raw(brokenLine)=[];
0054 
0055 line_key = regexprep(line_raw,'^ *-? ([^:]+)(:)($| .*)','$1');
0056 line_key = regexprep(line_key,'(.*!!omap)|(---)|( {4,}.*)','');
0057 
0058 line_value = regexprep(line_raw, '.*:$','');
0059 line_value = regexprep(line_value, '[^":]+: "?(.+)"?$','$1');
0060 line_value = regexprep(line_value, '(")|(^ {4,}- )','');
0061 line_value(strcmp(line_value,'''''')) = {''};
0062 line_value(strcmp(line_value,line_raw)) = {''};
0063 
0064 modelFields =   {'id',char();...
0065                'name',char();...
0066         'description',char();...
0067             'version',char();...
0068                'date',char();...
0069          'annotation',struct();...
0070                'rxns',{};...
0071            'rxnNames',{};...
0072                'mets',{};...
0073            'metNames',{};...
0074                   'S',sparse([]);...
0075                  'lb',{};... %Changed to double in the end.
0076                  'ub',{};... %Changed to double in the end.
0077                 'rev',{};... %Changed to double in the end.
0078                   'c',[];...
0079                   'b',cell(0,0);... %Changed to double in the end.
0080               'genes',cell(0,0);...
0081             'grRules',cell(0,0);...
0082          'rxnGeneMat',sparse([]);...
0083            'rxnComps',cell(0,0);... %Changed to double in the end.
0084          'subSystems',cell(0,0);...
0085             'eccodes',cell(0,0);...
0086          'rxnMiriams',cell(0,0);...
0087           'rxnDeltaG',{};... %Changed to double in the end.
0088            'rxnNotes',cell(0,0);...
0089       'rxnReferences',cell(0,0);...
0090 'rxnConfidenceScores',cell(0,0);...
0091            'metComps',cell(0,0);... %Changed to double in the end.
0092              'inchis',cell(0,0);...
0093           'metSmiles',cell(0,0);...
0094         'metFormulas',cell(0,0);...
0095          'metMiriams',cell(0,0);...
0096           'metDeltaG',{};... %Changed to double in the end.
0097          'metCharges',cell(0,0);... %Changed to double in the end.
0098            'metNotes',cell(0,0);...
0099               'comps',cell(0,0);...
0100           'compNames',cell(0,0);...
0101         'compOutside',cell(0,0);...
0102           'geneComps',cell(0,0);... %Changed to double in the end.
0103         'geneMiriams',cell(0,0);...
0104      'geneShortNames',cell(0,0);...
0105        'proteins',cell(0,0);...
0106       'unconstrained',cell(0,0);... %Changed to double in the end.
0107             'metFrom',cell(0,0);...
0108             'rxnFrom',cell(0,0)};
0109 for i=1:size(modelFields,1)
0110     model.(modelFields{i,1})=modelFields{i,2};
0111 end
0112 
0113 % If GECKO model
0114 if any(contains(line_key,'geckoLight'))
0115     isGECKO=true;
0116     ecFields = {'geckoLight', false;...
0117                       'rxns', {};...
0118                       'kcat', {};...
0119                     'source', cell(0,0);...
0120                      'notes', cell(0,0);...
0121                    'eccodes', cell(0,0);...
0122                      'genes', cell(0,0);...
0123                    'enzymes', cell(0,0);...
0124                         'mw', cell(0,0);...
0125                   'sequence', cell(0,0);...
0126                      'concs', cell(0,0);...
0127                  'rxnEnzMat', []};
0128     for i=1:size(ecFields,1)
0129         model.ec.(ecFields{i,1})=ecFields{i,2};
0130     end
0131     ecGecko=cell(25000,2);      ecGeckoNo=1;
0132     enzStoich=cell(100000,3);   enzStoichNo=1;
0133 else
0134     isGECKO=false;
0135 end
0136 
0137 section = 0;
0138 metMiriams=cell(100000,3);   metMirNo=1;
0139 rxnMiriams=cell(100000,3);   rxnMirNo=1;
0140 geneMiriams=cell(100000,3);  genMirNo=1;
0141 subSystems=cell(100000,2);   subSysNo=1;
0142 eccodes=cell(100000,2);      ecCodeNo=1;
0143 equations=cell(100000,3);    equatiNo=1;
0144 
0145 for i=1:numel(line_key)
0146     tline_raw = line_raw{i};
0147     tline_key = line_key{i};
0148     tline_value = line_value{i};
0149     % import different sections
0150     switch tline_raw
0151         case {'- metaData:', '- metaData: !!omap'}
0152             section = 1;
0153             if verbose
0154                 fprintf('\t%d\n', section);
0155             end
0156             continue % Go to next line
0157         case '- metabolites:'
0158             section = 2;
0159             if verbose
0160                 fprintf('\t%d\n', section);
0161             end
0162             pos=0;
0163             continue
0164         case '- reactions:'
0165             section = 3;
0166             if verbose
0167                 fprintf('\t%d\n', section);
0168             end
0169             pos=0;
0170             continue
0171         case {'- genes:', '- genes: []'}
0172             section = 4;
0173             if verbose
0174                 fprintf('\t%d\n', section);
0175             end
0176             pos=0;
0177             continue
0178         case '- compartments: !!omap'
0179             section = 5;
0180             if verbose
0181                 fprintf('\t%d\n', section);
0182             end
0183             pos=0;
0184             continue
0185         case '- ec-rxns:'
0186             section = 6;
0187             if verbose
0188                 fprintf('\t%d\n', section);
0189             end
0190             pos=0;
0191             continue
0192         case '- ec-enzymes:'
0193             section = 7;
0194             if verbose
0195                 fprintf('\t%d\n', section);
0196             end
0197             pos=0;
0198             continue
0199     end
0200 
0201     % skip over empty keys
0202     if isempty(tline_raw) || (isempty(tline_key) && contains(tline_raw,'!!omap'))
0203         continue;
0204     end
0205     
0206     % import metaData
0207     if section == 1
0208         switch tline_key
0209             case {'short_name','id'} %short_name used by human-GEM
0210                 model.id = tline_value;
0211             case 'name'
0212                 model.name = tline_value;                
0213             case 'full_name' %used by human-GEM
0214                 model.description = tline_value;
0215             case 'version'
0216                 model.version = tline_value;
0217             case 'date'
0218                 model.date = tline_value;                
0219             case 'taxonomy'
0220                 model.annotation.taxonomy = tline_value;
0221             case {'description','note'} %description used by human-GEM
0222                 model.annotation.note = tline_value;
0223             case 'github'
0224                 model.annotation.sourceUrl = tline_value;
0225             case 'sourceUrl'
0226                 model.annotation.sourceUrl = tline_value;                
0227             case 'givenName'
0228                 model.annotation.givenName = tline_value;
0229             case 'familyName'
0230                 model.annotation.familyName = tline_value;
0231             case 'authors'
0232                 model.annotation.authorList = tline_value;
0233             case 'email'
0234                 model.annotation.email = tline_value;
0235             case 'organization'
0236                 model.annotation.organization = tline_value;
0237             case {'geckoLight', 'gecko_light'}
0238                 if strcmp(tline_value,'true')
0239                     model.ec.geckoLight = true;
0240                 end
0241         end; continue
0242     end
0243 
0244     % import metabolites:
0245     if section == 2
0246         switch tline_key
0247             case 'id'
0248                 pos = pos + 1;
0249                 model = readFieldValue(model, 'mets', tline_value,pos);
0250                 readList=''; miriamKey='';
0251             case 'name'
0252                 model = readFieldValue(model, 'metNames', tline_value, pos);
0253                 readList=''; miriamKey='';
0254             case 'compartment'
0255                 model = readFieldValue(model, 'metComps', tline_value, pos);
0256                 readList=''; miriamKey='';
0257             case 'formula'
0258                 model = readFieldValue(model, 'metFormulas', tline_value, pos);
0259                 readList=''; miriamKey='';
0260             case 'charge'
0261                 model = readFieldValue(model, 'metCharges', tline_value, pos);
0262                 readList=''; miriamKey='';
0263             case 'notes'
0264                 model = readFieldValue(model, 'metNotes', tline_value, pos);
0265                 readList=''; miriamKey='';
0266             case 'inchis'
0267                 model = readFieldValue(model, 'inchis', tline_value, pos);
0268                 readList=''; miriamKey='';
0269             case 'smiles'
0270                 if isempty(tline_value)
0271                     readList = 'smilesInAnnotation';
0272                 else
0273                     model = readFieldValue(model, 'metSmiles', tline_value, pos);
0274                     readList=''; miriamKey='';
0275                 end    
0276             case 'deltaG'
0277                 model = readFieldValue(model, 'metDeltaG', tline_value, pos);
0278                 readList=''; miriamKey='';                                
0279             case 'metFrom'
0280                 model = readFieldValue(model, 'metFrom', tline_value, pos);
0281                 readList=''; miriamKey='';
0282             case 'annotation'
0283                 readList = 'annotation';
0284             otherwise
0285                 switch readList
0286                     case 'smilesInAnnotation'
0287                         %As for ec-code on a reaction: a bare list item
0288                         %continues the list, anything else ends it.
0289                         if isempty(tline_key)
0290                             model = readFieldValue(model, 'metSmiles', ...
0291                                 regexprep(tline_value,'^ +- "?(.*)"?$','$1'), pos);
0292                         else
0293                             readList = 'annotation';
0294                             [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo);
0295                         end
0296                     case 'annotation'
0297                         [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo);
0298                     otherwise
0299                         error(['Unknown entry in yaml file: ' tline_raw])
0300                 end
0301         end; continue
0302     end
0303 
0304     % import reactions:
0305     if section == 3
0306         switch tline_key
0307             case 'id'
0308                 pos = pos + 1;
0309                 model = readFieldValue(model, 'rxns', tline_value,pos);
0310                 readList=''; miriamKey='';
0311             case 'name'
0312                 model = readFieldValue(model, 'rxnNames', tline_value, pos);
0313                 readList=''; miriamKey='';
0314             case 'lower_bound'
0315                 model.lb(pos,1) = {tline_value};
0316                 readList=''; miriamKey='';
0317             case 'upper_bound'
0318                 model.ub(pos,1) = {tline_value};
0319                 readList=''; miriamKey='';
0320             case 'rev'
0321                 model.rev(pos,1) = {tline_value};
0322                 readList=''; miriamKey='';
0323             case 'gene_reaction_rule'
0324                 model = readFieldValue(model, 'grRules', tline_value, pos);
0325                 readList=''; miriamKey='';
0326             case {'rxnNotes', 'notes'}
0327                 model = readFieldValue(model, 'rxnNotes', tline_value, pos);
0328                 readList=''; miriamKey='';
0329             case 'rxnFrom'
0330                 model = readFieldValue(model, 'rxnFrom', tline_value, pos);
0331                 readList=''; miriamKey='';
0332             case 'deltaG'
0333                 model = readFieldValue(model, 'rxnDeltaG', tline_value, pos);
0334                 readList=''; miriamKey='';                
0335             case 'objective_coefficient'
0336                 model.c(pos,1) = 1;
0337                 readList=''; miriamKey='';
0338             case 'references'
0339                 model = readFieldValue(model, 'rxnReferences', tline_value, pos);
0340                 readList=''; miriamKey='';
0341             case 'confidence_score'
0342                 model = readFieldValue(model, 'rxnConfidenceScores', tline_value, pos);
0343                 readList=''; miriamKey='';
0344             case 'eccodes'
0345                 if isempty(tline_value)
0346                     readList = 'eccodes';
0347                 else
0348                     eccodes(ecCodeNo,1:2)={pos,tline_value};
0349                     ecCodeNo=ecCodeNo+1;
0350                 end
0351             case 'subsystem'
0352                 if isempty(tline_value)
0353                     readList = 'subsystem';
0354                 else
0355                     subSystems(subSysNo,1:2)={pos,tline_value};
0356                     subSysNo=subSysNo+1;
0357                 end
0358             case 'metabolites'
0359                 readList = 'equation';
0360             case 'ec-code'
0361                 %EC numbers emitted inside the annotation block rather than
0362                 %under the legacy top-level 'eccodes' key, which is still
0363                 %accepted above. A block list follows when the value is
0364                 %empty; it is read under its own state so that gathering can
0365                 %resume for the rest of the annotation block once the list
0366                 %ends, instead of swallowing every later entry as an EC code.
0367                 if isempty(tline_value)
0368                     readList = 'eccodesInAnnotation';
0369                 else
0370                     eccodes(ecCodeNo,1:2)={pos,tline_value};
0371                     ecCodeNo=ecCodeNo+1;
0372                 end
0373             case 'annotation'
0374                 readList = 'annotation';
0375                 
0376             otherwise
0377                 switch readList
0378                     case 'eccodes'
0379                         eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
0380                         ecCodeNo=ecCodeNo+1;
0381                     case 'eccodesInAnnotation'
0382                         %A bare list item ("- 1.1.2.4") has no colon, so
0383                         %tline_key is empty and the list continues. Anything
0384                         %else is the annotation block's next entry, so the
0385                         %list has ended.
0386                         if isempty(tline_key)
0387                             eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
0388                             ecCodeNo=ecCodeNo+1;
0389                         else
0390                             readList = 'annotation';
0391                             [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo);
0392                         end
0393                     case 'subsystem'
0394                         subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
0395                         subSysNo=subSysNo+1;
0396                     case 'annotation'
0397                         [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo);
0398                     case 'equation'
0399                         coeff = sscanf(tline_value,'%f');
0400                         equations(equatiNo,1:3)={pos,tline_key,coeff};
0401                         equatiNo=equatiNo+1;
0402                     otherwise
0403                         error(['Unknown entry in yaml file: ' tline_raw])
0404                 end                
0405         end; continue
0406     end
0407 
0408     % import genes:
0409     if section == 4
0410         switch tline_key
0411             case 'id'
0412                 pos = pos + 1;
0413                 model = readFieldValue(model, 'genes', tline_value, pos);
0414                 readList = '';
0415                 miriamKey = '';
0416             case 'name'
0417                 model = readFieldValue(model, 'geneShortNames', tline_value, pos);
0418             case 'protein'
0419                 model = readFieldValue(model, 'proteins', tline_value, pos);
0420             case 'annotation'
0421                 readList = 'annotation';
0422             otherwise
0423                 switch readList
0424                     case 'annotation'
0425                         [geneMiriams, miriamKey,genMirNo] = gatherAnnotation(pos,geneMiriams,tline_key,tline_value,miriamKey,genMirNo);
0426                     otherwise
0427                         error(['Unknown entry in yaml file: ' tline_raw])
0428                 end
0429         end; continue
0430     end
0431 
0432     % import compartments:
0433     if section == 5
0434         model.comps(end+1,1) = {tline_key};
0435         model.compNames(end+1,1) = {tline_value};
0436     end
0437 
0438     % import ec reaction info
0439     if section == 6
0440         switch tline_key
0441             case 'id'
0442                 pos = pos + 1;
0443                 model.ec = readFieldValue(model.ec, 'rxns', tline_value, pos);
0444                 readList='';
0445             case 'kcat'
0446                 model.ec = readFieldValue(model.ec, 'kcat', tline_value, pos);
0447                 readList='';
0448             case 'source'
0449                 model.ec = readFieldValue(model.ec, 'source', tline_value, pos);
0450                 readList='';
0451             case 'notes'
0452                 model.ec = readFieldValue(model.ec, 'notes', tline_value, pos);
0453                 readList='';
0454             case 'eccodes'
0455                 if isempty(tline_value)
0456                     readList = 'eccodes';
0457                 else
0458                     ecGecko(ecGeckoNo,1:2)={pos,tline_value};
0459                     ecGeckoNo=ecGeckoNo+1;
0460                 end
0461             case 'enzymes'
0462                 readList = 'enzStoich';
0463             otherwise
0464                 switch readList
0465                     case 'eccodes'
0466                         ecGecko(ecGeckoNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
0467                         ecGeckoNo=ecGeckoNo+1;
0468                     case 'enzStoich'
0469                         coeff = sscanf(tline_value,'%f');
0470                         enzStoich(enzStoichNo,1:3)={pos,tline_key,coeff};
0471                         enzStoichNo=enzStoichNo+1;
0472                     otherwise
0473                         error(['Unknown entry in yaml file: ' tline_raw])
0474                 end
0475         end; continue
0476     end
0477 
0478     % import ec enzyme info
0479     if section == 7
0480         switch tline_key
0481             case 'genes'
0482                 pos = pos + 1;
0483                 model.ec = readFieldValue(model.ec, 'genes', tline_value, pos);
0484             case 'enzymes'
0485                 model.ec = readFieldValue(model.ec, 'enzymes', tline_value, pos);
0486             case 'mw'
0487                 model.ec = readFieldValue(model.ec, 'mw', tline_value, pos);
0488             case 'sequence'
0489                 model.ec = readFieldValue(model.ec, 'sequence', tline_value, pos);
0490             case 'concs'
0491                 model.ec = readFieldValue(model.ec, 'concs', tline_value, pos);
0492             otherwise
0493                 error(['Unknown entry in yaml file: ' tline_raw])
0494         end; continue
0495     end
0496 end
0497 
0498 %Parse annotations
0499 if ~isempty(metMiriams)
0500     locs = cell2mat(metMiriams(:,1));
0501     for i=unique(locs)'
0502         model.metMiriams{i,1}.name=metMiriams(locs==i,2);
0503         model.metMiriams{i,1}.value=metMiriams(locs==i,3);
0504     end
0505 end
0506 if ~isempty(rxnMiriams)
0507     locs = cell2mat(rxnMiriams(:,1));
0508     for i=unique(locs)'
0509         model.rxnMiriams{i,1}.name=rxnMiriams(locs==i,2);
0510         model.rxnMiriams{i,1}.value=rxnMiriams(locs==i,3);
0511     end
0512 end
0513 if ~isempty(geneMiriams)
0514     locs = cell2mat(geneMiriams(:,1));
0515     for i=unique(locs)'
0516         model.geneMiriams{i,1}.name=geneMiriams(locs==i,2);
0517         model.geneMiriams{i,1}.value=geneMiriams(locs==i,3);
0518     end
0519 end
0520 
0521 %Parse subSystems
0522 if ~isempty(subSystems)
0523     locs = cell2mat(subSystems(:,1));
0524     for i=unique(locs)'
0525         model.subSystems{i,1}=subSystems(locs==i,2);
0526     end
0527 end
0528 
0529 %Parse ec-codes
0530 if ~isempty(eccodes)
0531     locs = cell2mat(eccodes(:,1));
0532     for i=unique(locs)'
0533         eccodesCat=strjoin(eccodes(locs==i,2),';');
0534         model.eccodes{i,1}=eccodesCat;
0535     end
0536     emptyEc=cellfun('isempty',model.eccodes);
0537     model.eccodes(emptyEc)={''};
0538 end
0539 
0540 % follow-up data processing
0541 if verbose
0542     fprintf('\nimporting completed\nfollow-up processing...');
0543 end
0544 [~, model.metComps] = ismember(model.metComps, model.comps);
0545 [~, model.geneComps] = ismember(model.geneComps, model.comps);
0546 [~, model.rxnComps] = ismember(model.rxnComps, model.comps);
0547 
0548 % Fill S-matrix
0549 rxnIdx = cellfun('isempty', equations(:,1));
0550 equations(rxnIdx,:) = '';
0551 rxnIdx = cell2mat(equations(:,1));
0552 [~,metIdx] = ismember(equations(:,2),model.mets);
0553 coeffs = cell2mat(equations(:,3));
0554 model.S=sparse(max(metIdx),max(rxnIdx));
0555 linearIndices = sub2ind([max(metIdx), max(rxnIdx)],metIdx,rxnIdx);
0556 model.S(linearIndices) = coeffs;
0557 
0558 % Convert strings to numeric
0559 model.metCharges = str2double(model.metCharges);
0560 model.lb = str2double(model.lb);
0561 model.ub = str2double(model.ub);
0562 model.rxnConfidenceScores = str2double(model.rxnConfidenceScores);
0563 model.b = zeros(length(model.mets),1);
0564 model.metDeltaG = str2double(model.metDeltaG);
0565 model.rxnDeltaG = str2double(model.rxnDeltaG);
0566 
0567 % Fill some other fields
0568 model.annotation.defaultLB = min(model.lb);
0569 model.annotation.defaultUB = max(model.ub);
0570 if numel(model.lb)<numel(model.rxns) %No LB reported = min
0571     model.lb(end+1:numel(model.rxns)-numel(model.lb),1) = double(model.annotation.defaultLB);
0572 end
0573 if numel(model.ub)<numel(model.rxns) %No UB reported = max
0574     model.ub(end+1:numel(model.rxns)-numel(model.ub),1) = double(model.annotation.defaultUB);
0575 end
0576 if ~all(cellfun('isempty',model.rev))
0577     model.rev = str2double(model.rev);
0578 else
0579     model.rev = [];
0580 end
0581 if numel(model.rev)<numel(model.rxns) %No rev reported, assume from LB and UB
0582     model.rev(end+1:numel(model.rxns)-numel(model.rev),1) = double(model.lb<0 & model.ub>0);
0583 end
0584 
0585 % Remove empty fields, otherwise fill to correct length
0586 % Reactions
0587 for i={'rxnNames','grRules','eccodes','rxnNotes','rxnReferences',...
0588        'rxnFrom','subSystems','rxnMiriams'} % Empty strings
0589    model = emptyOrFill(model,i{1},{''},'rxns');
0590 end
0591 for i={'c'} % Zeros
0592    model = emptyOrFill(model,i{1},0,'rxns',true);
0593 end
0594 for i={'rxnConfidenceScores','rxnDeltaG'} % NaNs
0595    model = emptyOrFill(model,i{1},NaN,'rxns');
0596 end
0597 for i={'rxnComps'} % Ones, assume first compartment
0598    model = emptyOrFill(model,i{1},1,'rxns');
0599 end
0600 % Metabolites
0601 for i={'metNames','inchis','metFormulas','metMiriams','metFrom','metSmiles','metNotes'} % Empty strings
0602    model = emptyOrFill(model,i{1},{''},'mets');
0603 end
0604 for i={'metCharges','unconstrained'} % Zeros
0605    model = emptyOrFill(model,i{1},0,'mets');
0606 end
0607 for i={'metDeltaG'} % % NaNs
0608     model = emptyOrFill(model,i{1},NaN,'mets');
0609  end
0610 for i={'metComps'} % Ones, assume first compartment
0611    model = emptyOrFill(model,i{1},1,'mets');
0612 end
0613 % Genes
0614 for i={'geneMiriams','geneShortNames','proteins'} % Empty strings
0615    model = emptyOrFill(model,i{1},{''},'genes');
0616 end
0617 for i={'geneComps'} % Ones, assume first compartment
0618    model = emptyOrFill(model,i{1},1,'genes');
0619 end
0620 % Comps
0621 for i={'compNames'} % Empty strings
0622    model = emptyOrFill(model,i{1},{''},'comps');
0623 end
0624 for i={'compOutside'} % First comp
0625    model = emptyOrFill(model,i{1},model.comps{1},'comps');
0626 end
0627 % Single fields are kept, even if empty
0628 % for i={'description','name','version','date','annotation'}
0629 %     if isempty(model.(i{1}))
0630 %         model = rmfield(model,i{1});
0631 %     end
0632 % end
0633 
0634 if isfield(model,'subSystems')
0635 % If all entries are 1x1, then flatten
0636     if all(cellfun(@(x) numel(x) <= 1, model.subSystems))
0637         model.subSystems = transpose([model.subSystems{:}]);
0638     end
0639 end
0640 
0641 % Make rxnGeneMat fields and map to the existing model.genes field
0642 if isfield(model,'grRules')
0643     [genes, rxnGeneMat] = getGenesFromGrRules(model.grRules);
0644     model.rxnGeneMat = sparse(numel(model.rxns),numel(model.genes));
0645     [~,geneOrder] = ismember(genes,model.genes);
0646     if any(geneOrder == 0)
0647         error(['The grRules includes the following gene(s), that are not in '...
0648             'the list of model genes: ', genes{~geneOrder}])
0649     end
0650     model.rxnGeneMat(:,geneOrder) = rxnGeneMat;
0651 end
0652 
0653 % Finalize GECKO model
0654 if isGECKO
0655     % Fill in empty fields and empty entries
0656     for i={'kcat','source','notes','eccodes'} % Even keep empty
0657         model.ec = emptyOrFill(model.ec,i{1},{''},'rxns',true);
0658     end
0659     for i={'enzymes','mw','sequence'}
0660         model.ec = emptyOrFill(model.ec,i{1},{''},'genes',true);
0661     end
0662     model.ec = emptyOrFill(model.ec,'concs',{'NaN'},'genes',true);
0663     model.ec = emptyOrFill(model.ec,'kcat',{'0'},'genes',true);
0664     % Change string to double
0665     for i={'kcat','mw','concs'}
0666         if isfield(model.ec,i{1})
0667             model.ec.(i{1}) = str2double(model.ec.(i{1}));
0668         end
0669     end
0670     % Fill rxnEnzMat
0671     rxnIdx              = cellfun('isempty', enzStoich(:,1));
0672     enzStoich(rxnIdx,:) = '';
0673     rxnIdx              = cell2mat(enzStoich(:,1));
0674     [~,enzIdx]          = ismember(enzStoich(:,2),model.ec.enzymes);
0675     coeffs              = cell2mat(enzStoich(:,3));
0676     model.ec.rxnEnzMat  = zeros(numel(model.ec.rxns), numel(model.ec.genes));
0677     linearIndices       = sub2ind([numel(model.ec.rxns), numel(model.ec.genes)], rxnIdx, enzIdx);
0678     model.ec.rxnEnzMat(linearIndices) = coeffs;
0679     %Parse ec-codes
0680     if ~isempty(ecGecko)
0681         locs = cell2mat(ecGecko(:,1));
0682         for i=unique(locs)'
0683             ecGeckoCat=strjoin(ecGecko(locs==i,2),';');
0684             model.ec.eccodes{i,1}=ecGeckoCat;
0685         end
0686         emptyEc=cellfun('isempty',model.ec.eccodes);
0687         model.ec.eccodes(emptyEc)={''};
0688     end
0689 end
0690 
0691 if verbose
0692     fprintf(' Done!\n');
0693 end
0694 end
0695 
0696 function model = emptyOrFill(model,field,emptyEntry,type,keepEmpty)
0697 if nargin<5
0698     keepEmpty=false;
0699 end
0700 if isnumeric(emptyEntry)
0701     emptyCells=isempty(model.(field));
0702 else
0703     emptyCells=cellfun('isempty',model.(field));
0704 end
0705 if all(emptyCells) && ~keepEmpty
0706     model = rmfield(model, field);
0707 elseif numel(model.(field))<numel(model.(type))
0708     model.(field)(end+1:numel(model.(type)),1)=emptyEntry;
0709 end
0710 end
0711 
0712 function model = readFieldValue(model, fieldName, value, pos)
0713 if numel(model.(fieldName))<pos-1
0714     model.(fieldName)(end+1:pos,1) = {''};
0715 end
0716 model.(fieldName)(pos,1) = {value};
0717 end
0718 
0719 function [miriams,miriamKey,entryNumber] = gatherAnnotation(pos,miriams,key,value,miriamKey,entryNumber)
0720 if isempty(key)
0721     key=miriamKey;
0722 else
0723     miriamKey=key;
0724 end
0725 if ~isempty(value)
0726     miriams(entryNumber,1:3) = {pos, key, strip(value)};
0727     entryNumber = entryNumber + 1;
0728 end
0729 end