fix: read RAVEN 3 model files, and tolerate a non-double rev

INIT/ftINIT.m

@@ -135,7 +135,7 @@
         %Get the matching mets
         metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i}));
         prodRxnsSel = any(prepData.refModel.S(metSel,:) > 0,1) | ... %direct producers
-                     (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.'); %reversible reactions that are consumers
+                     (any(prepData.refModel.S(metSel,:) < 0,1) & prepData.refModel.rev.' ~= 0); %reversible reactions that are consumers
         %convert the production rxns from refModel to minModel
         prepData.groupIds;
         [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns);

INIT/mergeLinear.m

@@ -62,8 +62,12 @@
     for i=1:numel(twoNonZero)
         involvedRxns=find(reducedModel.S(twoNonZero(i),:));
         %Check that we can have one positive and one negative
-        pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns));
-        neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns));
+        %rev is compared rather than used as-is: the other operand is a
+        %sparse logical, and MATLAB refuses to combine one with an integer
+        %array ("Sparse integer array arithmetic operations are not
+        %supported"). A model built outside RAVEN can carry an integer rev.
+        pos = sum(reducedModel.S(twoNonZero(i),involvedRxns).' > 0 | reducedModel.rev(involvedRxns) ~= 0);
+        neg = sum(reducedModel.S(twoNonZero(i),involvedRxns).' < 0 | reducedModel.rev(involvedRxns) ~= 0);
 
 
         %Check so that one or both of the reactions haven't been merged

doc/INIT/ftINIT.html

@@ -270,7 +270,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
 0135         <span class="comment">%Get the matching mets</span>
 0136         metSel = ismember(upper(prepData.refModel.metNames),upper(metabolomicsData{i}));
 0137         prodRxnsSel = any(prepData.refModel.S(metSel,:) &gt; 0,1) | <span class="keyword">...</span><span class="comment"> %direct producers</span>
-0138                      (any(prepData.refModel.S(metSel,:) &lt; 0,1) &amp; prepData.refModel.rev.'); <span class="comment">%reversible reactions that are consumers</span>
+0138                      (any(prepData.refModel.S(metSel,:) &lt; 0,1) &amp; prepData.refModel.rev.' ~= 0); <span class="comment">%reversible reactions that are consumers</span>
 0139         <span class="comment">%convert the production rxns from refModel to minModel</span>
 0140         prepData.groupIds;
 0141         [~,ia,ib] = intersect(prepData.minModel.rxns,prepData.refModel.rxns);

io/importModel.m

@@ -739,11 +739,22 @@
     end
 end
 if isfield(modelSBML,'notes')
-    startString=strfind(modelSBML.notes,'xhtml">');
-    endString=strfind(modelSBML.notes,'</body>');
-    if any(startString) && any(endString)
-        model.annotation.note=modelSBML.notes(startString(1)+7:endString-1);
-        model.annotation.note=regexprep(model.annotation.note,'<p.*?>|</p.*?>','');
+    %Two shapes reach here. RAVEN's own exportModel writes the note as bare
+    %text inside a <body>; a writer that serialises a notes dictionary emits
+    %one '<p>key: value</p>' per entry and no <body> at all, carrying the note
+    %under its own key. Take the labelled entry when there is one, so the
+    %other entries (version, and whatever else the dictionary held) stay out
+    %of the note, and fall back to the <body> text otherwise.
+    note=parseNote(modelSBML.notes,'note');
+    if isempty(note)
+        startString=strfind(modelSBML.notes,'xhtml">');
+        endString=strfind(modelSBML.notes,'</body>');
+        if any(startString) && any(endString)
+            note=modelSBML.notes(startString(1)+7:endString-1);
+        end
+    end
+    if ~isempty(note)
+        model.annotation.note=regexprep(note,'<p.*?>|</p.*?>','');
         model.annotation.note=strtrim(model.annotation.note);
         if regexp(model.annotation.note,'This file was generated using the exportModel function in RAVEN Toolbox \d\.\d and OutputSBML in libSBML')
             model.annotation=rmfield(model.annotation,'note'); % Default note added when running exportModel

io/readYAMLmodel.m

@@ -148,7 +148,7 @@
     tline_value = line_value{i};
     % import different sections
     switch tline_raw
-        case '- metaData:'
+        case {'- metaData:', '- metaData: !!omap'}
             section = 1;
             if verbose
                 fprintf('\t%d\n', section);
@@ -168,7 +168,7 @@
             end
             pos=0;
             continue
-        case '- genes:'
+        case {'- genes:', '- genes: []'}
             section = 4;
             if verbose
                 fprintf('\t%d\n', section);
@@ -234,7 +234,7 @@
                 model.annotation.email = tline_value;
             case 'organization'
                 model.annotation.organization = tline_value;
-            case 'geckoLight'
+            case {'geckoLight', 'gecko_light'}
                 if strcmp(tline_value,'true')
                     model.ec.geckoLight = true;
                 end
@@ -267,8 +267,12 @@
                 model = readFieldValue(model, 'inchis', tline_value, pos);
                 readList=''; miriamKey='';
             case 'smiles'
-                model = readFieldValue(model, 'metSmiles', tline_value, pos);
-                readList=''; miriamKey='';    
+                if isempty(tline_value)
+                    readList = 'smilesInAnnotation';
+                else
+                    model = readFieldValue(model, 'metSmiles', tline_value, pos);
+                    readList=''; miriamKey='';
+                end    
             case 'deltaG'
                 model = readFieldValue(model, 'metDeltaG', tline_value, pos);
                 readList=''; miriamKey='';                                
@@ -279,6 +283,16 @@
                 readList = 'annotation';
             otherwise
                 switch readList
+                    case 'smilesInAnnotation'
+                        %As for ec-code on a reaction: a bare list item
+                        %continues the list, anything else ends it.
+                        if isempty(tline_key)
+                            model = readFieldValue(model, 'metSmiles', ...
+                                regexprep(tline_value,'^ +- "?(.*)"?$','$1'), pos);
+                        else
+                            readList = 'annotation';
+                            [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo);
+                        end
                     case 'annotation'
                         [metMiriams, miriamKey, metMirNo] = gatherAnnotation(pos,metMiriams,tline_key,tline_value,miriamKey,metMirNo);
                     otherwise
@@ -309,7 +323,7 @@
             case 'gene_reaction_rule'
                 model = readFieldValue(model, 'grRules', tline_value, pos);
                 readList=''; miriamKey='';
-            case 'rxnNotes'
+            case {'rxnNotes', 'notes'}
                 model = readFieldValue(model, 'rxnNotes', tline_value, pos);
                 readList=''; miriamKey='';
             case 'rxnFrom'
@@ -343,6 +357,19 @@
                 end
             case 'metabolites'
                 readList = 'equation';
+            case 'ec-code'
+                %EC numbers emitted inside the annotation block rather than
+                %under the legacy top-level 'eccodes' key, which is still
+                %accepted above. A block list follows when the value is
+                %empty; it is read under its own state so that gathering can
+                %resume for the rest of the annotation block once the list
+                %ends, instead of swallowing every later entry as an EC code.
+                if isempty(tline_value)
+                    readList = 'eccodesInAnnotation';
+                else
+                    eccodes(ecCodeNo,1:2)={pos,tline_value};
+                    ecCodeNo=ecCodeNo+1;
+                end
             case 'annotation'
                 readList = 'annotation';
 
@@ -351,6 +378,18 @@
                     case 'eccodes'
                         eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
                         ecCodeNo=ecCodeNo+1;
+                    case 'eccodesInAnnotation'
+                        %A bare list item ("- 1.1.2.4") has no colon, so
+                        %tline_key is empty and the list continues. Anything
+                        %else is the annotation block's next entry, so the
+                        %list has ended.
+                        if isempty(tline_key)
+                            eccodes(ecCodeNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
+                            ecCodeNo=ecCodeNo+1;
+                        else
+                            readList = 'annotation';
+                            [rxnMiriams, miriamKey,rxnMirNo] = gatherAnnotation(pos,rxnMiriams,tline_key,tline_value,miriamKey,rxnMirNo);
+                        end
                     case 'subsystem'
                         subSystems(subSysNo,1:2)={pos,regexprep(tline_value,'^ +- "?(.*)"?$','$1')};
                         subSysNo=subSysNo+1;

readme/RAVEN_structure_fields.csv

@@ -49,7 +49,7 @@ geneMiriams;;Struct;Structure with MIRIAM information about the genes;Non-mandat
 ;value;Column Cell Array of Strings;The ID for particular MIRIAM entry;Non-mandatory field;P12345;
 geneShortNames;;Column Cell Array of Strings;Gene alternative names (e.g. ERG10);Non-mandatory field;ERG10;
 metNames;;Column Cell Array of Strings;Metabolite names;Non-mandatory field;alpha-D-glucose;
-metComps;;Column Vector of Doubles;Compartments for metabolites;Mandatory field;1;Shows the index of affiliated compartment in comps, not the compartment abbreviation
+metComps;;Column Vector of Integers;Compartments for metabolites;Mandatory field;1;Written as an integer so a reader that converts the 1-based index into a position gets one; MATLAB stores an integer-valued double narrowed to an integer type anyway, and both classes index comps identically.
 inchis;;Column Cell Array of Strings;InChI-codes for metabolites;Non-mandatory field;1/C6H12O6/c7-1-2-3(8)4(9)5(10)6(11)12-2/h2-11H,1H2/t2-,3-,4+,5-,6+/m1/s1;
 metFormulas;;Column Cell Array of Strings;Chemical formulas for metabolites;Non-mandatory field;C6H12O6;
 metMiriams;;Struct;Structure with MIRIAM information about the metabolites;Non-mandatory field;;If there are multiple miriams, 'name' indexes must match 'value' indexes