RAVEN 2.11.4

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

core/replaceMets.m

@@ -20,7 +20,7 @@
 % Output:
 %   model           model structure with selected metabolites replaced
 %   removedRxns     identifiers of duplicate reactions that were removed
-%   idxDuplRxns     index of removedRxns in original model
+%   idxDuplRxns     indices of removedRxns in the input model
 %
 % Note: This function is useful when the model contains both 'oxygen' and
 % 'o2' as metabolite names. If 'oxygen' and 'o2' are identifiers instead,
@@ -59,10 +59,11 @@
     error('The to-be-replaced metabolite cannot be found in the model.');
 end
 
-rxnsWithMet = find(model.S(metIdx,:));
+[~, rxnsWithMet] = find(model.S(metIdx,:));
+rxnsWithMet = unique(rxnsWithMet);
 if verbose==true
-    fprintf('\n\nThe following reactions contain the to-be-replaced metabolite as reactant:\n')
-    fprintf(strjoin(model.rxns(rxnsWithMet),'\n'))
+    fprintf('\n\nThe following reactions involve the to-be-replaced metabolite:\n')
+    fprintf('%s\n',model.rxns{rxnsWithMet})
     fprintf('\n')
 end
 
@@ -88,8 +89,11 @@
 
 idxDelete=[];
 if identifiers
+    if any(ismember(metIdx,repIdx))
+        error('The to-be-replaced metabolite and the replacement metabolite are the same.');
+    end
     originalStoch = model.S(metIdx,rxnsWithMet);
-    model.S(repIdx,rxnsWithMet) = originalStoch;
+    model.S(repIdx,rxnsWithMet) = model.S(repIdx,rxnsWithMet) + originalStoch;
     model.S(metIdx,rxnsWithMet) = 0;
     idxDelete = metIdx;
 else
@@ -99,54 +103,38 @@
     % the to-be-replace metabolite ID deleted.
 
     % Build list of metaboliteName[compartment]
-    metCompsN =cellstr(num2str(model.metComps));
-    map = containers.Map(cellstr(num2str(transpose(1:length(model.comps)))),model.comps);
-    metCompsN = map.values(metCompsN);
-    metCompsN = strcat(lower(model.metNames),'[',metCompsN,']');
+    comps = model.comps(:);
+    metCompsN = strcat(lower(model.metNames(:)),'[',comps(model.metComps(:)),']');
 
-    for i = 1:length(repIdx)
-        metCompsNidx=find(strcmp(metCompsN(repIdx(i)), metCompsN));
-        if length(metCompsNidx)>1
-            for j = 2:length(metCompsNidx)
-                model.S(metCompsNidx(1),:) = model.S(metCompsNidx(1),:) + model.S(metCompsNidx(j),:);
-                idxDelete=[idxDelete; metCompsNidx(j)]; % Make list of metabolite IDs to delete
-            end
+    % Within each metaboliteName[compartment] group, the first replacement
+    % metabolite is kept and the other members are merged into it.
+    groups = unique(metCompsN(repIdx),'stable');
+    for i = 1:numel(groups)
+        metCompsNidx = find(strcmp(groups{i}, metCompsN));
+        keepIdx = metCompsNidx(find(ismember(metCompsNidx,repIdx),1));
+        others = metCompsNidx(metCompsNidx ~= keepIdx);
+        if ~isempty(others)
+            model.S(keepIdx,:) = model.S(keepIdx,:) + sum(model.S(others,:),1);
+            idxDelete = [idxDelete; others]; % Make list of metabolite IDs to delete
         end
     end
 end
-
 if ~isempty(idxDelete)
-    model.S(idxDelete,:) =[];
-    model.mets(idxDelete) = [];
-    model.metNames(idxDelete) = [];
-    model.metComps(idxDelete) = [];
-    model.b(idxDelete) = [];
-    if isfield(model,'metFormulas')
-        model.metFormulas(idxDelete) = [];
-    end
-    if isfield(model,'unconstrained')
-        model.unconstrained(idxDelete) = [];
-    end
-    if isfield(model,'metMiriams')
-        model.metMiriams(idxDelete) = [];
-    end
-    if isfield(model,'metCharges')
-        model.metCharges(idxDelete) = [];
-    end
-    if isfield(model,'metDeltaG')
-        model.metDeltaG(idxDelete) = [];
-    end
-    if isfield(model,'inchis')
-        model.inchis(idxDelete) = [];
-    end
-    if isfield(model,'metSmiles')
-        model.metSmiles(idxDelete) = [];
-    end
-    if isfield(model,'metFrom')
-        model.metFrom(idxDelete) = [];
-    end
+    repMets = model.mets(repIdx);
+    model = removeMets(model,idxDelete);
+    % Indices of the replacement metabolites shift when metabolites with a
+    % lower index are deleted.
+    [~, repIdx] = ismember(repMets,model.mets);
+    repIdx(repIdx==0) = [];
 end
 
 % This could now have created duplicate reactions. Contract model.
-model=contractModel(model,[],repIdx);
+if isempty(idxDelete)
+    removedRxns = cell(0,1);
+    idxDuplRxns = zeros(0,1);
+else
+    originalRxns = model.rxns;
+    [model, removedRxns] = contractModel(model,[],repIdx);
+    [~, idxDuplRxns] = ismember(removedRxns,originalRxns);
+end
 end

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);