fix

findRootDeadEndMets.m

@@ -0,0 +1,43 @@
+function [justConsumedMets, justProducedMets]=findRootDeadEndMets(model)
+  %   This script will find metabolites which are just produced or just 
+  %   consumed in the reactions. These metabolites cause
+  %   deadEnd reactions.
+  %
+  %   model                 a model structure
+  %   justConsumedMets      list of metabolites which are just consumed by
+  %   the model
+  %   justProducedMets      list of metabolites which are just produced by
+  %   the model
+  %
+  %   Rasool Saghaleyni, 2020-03-07
+
+  % Adding all boundry metabolites to the model
+  model = addBoundaryMets(model);
+  % remove all both directional reactions from model
+  uniModel=removeReactions(model,model.rxns(model.rev==1));
+  % generate a model just containing both directional reactions
+  forwModel=removeReactions(model,model.rxns(model.rev==0));
+  % fixing directionality of direction
+  forwModel.rev(:,1) = 0;
+  % generate a model just containing both directional reactions
+  revModel=removeReactions(model,model.rxns(model.rev==0));
+  % change directionality of reaction to reverse
+  revModel.S=-revModel.S;
+  % fixing directionality of direction
+  revModel.rev(:,1) = 0;
+  % adding removed reactions to main model
+  newModel=mergeModels({uniModel,forwModel,revModel});
+  % finding metabolites which are just produced or consumed in all
+  % reactions
+  for i=1:size(newModel.mets)
+    if size(find(newModel.S(i,:)<0),2)==0
+        all_mets(i,1)={'prod'};
+    elseif size(find(newModel.S(i,:)>0),2)==0
+        all_mets(i,1)={'cons'};
+    end
+  end
+  % if you add this metabolites to the model using addExchangeRxns function, there will
+  % not be any deadEnd reaction in the model.
+  justConsumedMets=newModel.mets(strcmp(all_mets,'cons'));
+  justProducedMets=newModel.mets(strcmp(all_mets,'prod'));
+end

genericHumanSec.m

@@ -1,4 +1,4 @@
-function [ihumanSec, secRxns]=genericHumanSec(ihuman, tempRxns)
+function [ihumanSec, secAddon]=genericHumanSec(ihuman, tempRxns)
   % humanSec
   %   This script generates and adds all secretory reactions to generic metabolic model of human (ihuman)
   %
@@ -343,9 +343,9 @@
   newMets = setdiff(metList,outModel.metNamesC); %list of new metabolites with their compartment.
   if isempty(newMets) == 0
     [newNames, newComps] = splitComp(newMets);
-    metsToAdd.mets = newMets;
     metsToAdd.compartments = newComps;
     metsToAdd.metNames = newNames;
+    metsToAdd.mets = strcat(metsToAdd.metNames,metsToAdd.compartments);
     secModel = addMets(secModel,metsToAdd);
   end
   %------------------------------------------------------------------------------
@@ -355,19 +355,22 @@
   rxnsToAdd.rxnNames = rxnNames(ismember(rxnNames,secModel.rxns) == 0);
   rxnsToAdd.grRules = GPRs(ismember(rxnNames,secModel.rxns) == 0);
   rxnsToAdd.subSystems(1:numel(rxnsToAdd.rxnNames),1) = {'protein secretion'};
-
-  com = Comps(ismember(rxnNames,secModel.rxns) == 0);
+  %------------------------------------------------------------------------------
   %Add protein secretion reactions to human1-drived model and modify grule Matrix
   secModel = addRxns(secModel,rxnsToAdd,3,[],false);
-  %Add compartments for new reactions
-  %secModel.rxnComps(numel(model.rxns)+1:numel(secModel.rxns)) = com(1:end);
-  ihumanSec = secModel;
   %------------------------------------------------------------------------------
-  addedRxns.rxnsNames = rxnsToAdd.rxns;
-  addedRxns.rxnFormula = rxnsToAdd.equations;
-  addedRxns.rxnGPRs = rxnsToAdd.grRules;
-  addedRxns.rxnComps = com;
-  addedRxns.newMets = newMets;
-  addedRxns.newGens = genesToAdd.genes;
-  secRxns = addedRxns;
+  [justConsumedMets, justProducedMets]=findRootDeadEndMets(secModel);
+  justConsumedMets_sec = intersect(metsToAdd.mets,justConsumedMets);
+  justProducedMets_sec = intersect(metsToAdd.mets,justProducedMets);
+  if size(justConsumedMets_sec,1)~=0
+    [secModel, ~ ]=addExchangeRxns(secModel, 'in', justConsumedMets_sec);
+  end
+  if size(justProducedMets_sec,1)~=0
+    [secModel, ~ ]=addExchangeRxns(secModel, 'out', justProducedMets_sec);
+  end
+  ihumanSec=secModel;
+  %------------------------------------------------------------------------------
+  rxns2rm=intersect(ihumanSec.rxns,ihuman.rxns);
+  secAddon=removeReactions(ihumanSec,rxns2rm,1,1,1);
+  %------------------------------------------------------------------------------
 end

humanSec.m

@@ -394,9 +394,9 @@
   newMets = setdiff(metList,outModel.metNamesC); %list of new metabolites with their compartment.
   if isempty(newMets) == 0
     [newNames, newComps] = splitComp(newMets);
-    metsToAdd.mets = newMets;
     metsToAdd.compartments = newComps;
     metsToAdd.metNames = newNames;
+    metsToAdd.mets = strcat(metsToAdd.metNames,metsToAdd.compartments);
     secModel = addMets(secModel,metsToAdd);
   end
   %------------------------------------------------------------------------------