addTransport

 addTransport
   Adds transport reactions between compartments

   model           a model structure
   fromComp        the id of the compartment to transport from (should
                   match model.comps)
   toComps         a cell array of compartment names to transport to (should
                   match model.comps)
   metNames        the metabolite names to add transport for (optional, all
                   metabolites in fromComp)
   isRev           true if the transport reactions should be reversible
                   (optional, default true)
   onlyToExisting  true if transport of a metabolite should only be added
                   if it already exists in toComp. If false, then new metabolites
                   are added with addMets first (optional, default true)
   prefix          string specifying prefix to reaction IDs (optional, default
                   'tr_')

   This is a faster version than addRxns when adding transport reactions.
   New reaction names are formatted as "metaboliteName, fromComp-toComp", 
   while new reaction IDs are sequentially counted with a tr_ prefix:
   e.g. tr_0001, tr_0002, etc.

 Usage: [model, addedRxns]=addTransport(model,fromComp,toComps,metNames,...
           isRev,onlyToExisting,prefix)

0001 function [model, addedRxns]=addTransport(model,fromComp,toComps,metNames,isRev,onlyToExisting,prefix)
0002 % addTransport
0003 %   Adds transport reactions between compartments
0004 %
0005 %   model           a model structure
0006 %   fromComp        the id of the compartment to transport from (should
0007 %                   match model.comps)
0008 %   toComps         a cell array of compartment names to transport to (should
0009 %                   match model.comps)
0010 %   metNames        the metabolite names to add transport for (optional, all
0011 %                   metabolites in fromComp)
0012 %   isRev           true if the transport reactions should be reversible
0013 %                   (optional, default true)
0014 %   onlyToExisting  true if transport of a metabolite should only be added
0015 %                   if it already exists in toComp. If false, then new metabolites
0016 %                   are added with addMets first (optional, default true)
0017 %   prefix          string specifying prefix to reaction IDs (optional, default
0018 %                   'tr_')
0019 %
0020 %   This is a faster version than addRxns when adding transport reactions.
0021 %   New reaction names are formatted as "metaboliteName, fromComp-toComp",
0022 %   while new reaction IDs are sequentially counted with a tr_ prefix:
0023 %   e.g. tr_0001, tr_0002, etc.
0024 %
0025 % Usage: [model, addedRxns]=addTransport(model,fromComp,toComps,metNames,...
0026 %           isRev,onlyToExisting,prefix)
0027 
0028 fromComp=char(fromComp);
0029 [I, fromID]=ismember(model.comps,fromComp);
0030 fromID=find(fromID);
0031 if sum(I)~=1
0032     EM='fromComps must have exactly one match in model.comps';
0033     dispEM(EM);
0034 end
0035 toComps=convertCharArray(toComps);
0036 [I, toIDs]=ismember(toComps,model.comps);
0037 if ~all(I)
0038     EM='All compartments in toComps must have a match in model.comps';
0039     dispEM(EM);
0040 end
0041 if nargin<4
0042     %Find all metabolites in fromComp
0043     metNames=model.metNames(model.metComps==fromID);
0044 elseif isempty(metNames)
0045     %Find all metabolites in fromComp
0046     metNames=model.metNames(ismember(model.metComps,model.comps(fromID)));
0047 else
0048     metNames=convertCharArray(metNames);
0049 end
0050 
0051 if nargin<5
0052     isRev=true;
0053 end
0054 if nargin<6
0055     onlyToExisting=true;
0056 end
0057 if nargin<7
0058     prefix='tr_';
0059 else
0060     prefix=char(prefix);
0061 end
0062 
0063 %Check that the names are unique
0064 if numel(unique(metNames))~=numel(metNames)
0065     dispEM('Not all metabolite names are unique');
0066 end
0067 
0068 %Get the indexes of the mets in fromComp
0069 I=find(model.metComps==fromID);
0070 [J, K]=ismember(metNames,model.metNames(I));
0071 if ~all(J)
0072     EM='Not all metabolites in metNames exist in fromComp';
0073     dispEM(EM);
0074 end
0075 fromMets=I(K); %These are the ids of the metabolites to transport. The order corresponds to metNames
0076 addedRxns={};
0077 %Loop through and add for each compartment in toComps
0078 for i=1:numel(toComps)
0079     fromMetsInComp=fromMets; %If onlyToExisting==true then not all mets are transported to each compartment
0080     %Get the indexes of the mets in the compartment
0081     I=find(model.metComps==toIDs(i));
0082     [J, K]=ismember(metNames,model.metNames(I));
0083     if onlyToExisting==true || all(J)
0084         toMets=I(K(J)); %Only look at the existing ones
0085         fromMetsInComp=fromMetsInComp(J);
0086     else
0087         %This is if not all metabolites exist in the target compartment,
0088         %and they should be added
0089         metsToAdd.metNames=metNames(J==0);
0090         metsToAdd.compartments=toComps{i};
0091         model=addMets(model,metsToAdd);
0092         
0093         %Redo the mapping when all mets are there. A bit lazy, but it's
0094         %fast anyways
0095         I=find(model.metComps==toIDs(i));
0096         [~, K]=ismember(metNames,model.metNames(I));
0097         toMets=I(K); %All are guaranteed to be found now
0098     end
0099     
0100     %Construct the S matrix
0101     nRxns=numel(fromMetsInComp);
0102     newS=zeros(numel(model.mets),nRxns);
0103     newS(sub2ind(size(newS),fromMetsInComp(:),(1:nRxns)'))=-1;
0104     newS(sub2ind(size(newS),toMets(:),(1:nRxns)'))=1;
0105     
0106     %Add the reactions
0107     model.S=[model.S sparse(newS)];
0108     if isfield(model.annotation,'defaultLB')
0109         lb = model.annotation.defaultLB;
0110         ub = model.annotation.defaultUB;
0111     else
0112         lb = -inf;
0113         ub = inf;
0114     end
0115     if isRev==true
0116         model.lb=[model.lb;ones(nRxns,1)*lb];
0117         model.rev=[model.rev;ones(nRxns,1)];
0118     else
0119         model.lb=[model.lb;zeros(nRxns,1)];
0120         model.rev=[model.rev;zeros(nRxns,1)];
0121     end
0122     model.ub=[model.ub;ones(nRxns,1)*ub];
0123     model.c=[model.c;zeros(nRxns,1)];
0124     
0125     %Add annotation
0126     filler=cell(nRxns,1);
0127     filler(:)={''};
0128     addedRxnsID=generateNewIds(model,'rxns',prefix,nRxns);
0129     addedRxnsName=transpose(strcat(metNames, {' transport, '}, model.compNames(fromID), '-', model.compNames(toIDs(i))));
0130     model.rxns=[model.rxns;addedRxnsID];
0131     model.rxnNames=[model.rxnNames;addedRxnsName];
0132     
0133     if isfield(model,'eccodes')
0134         model.eccodes=[model.eccodes;filler];
0135     end
0136     if isfield(model,'subSystems')
0137         ssFiller=filler;
0138         if isRev==1
0139             ssFiller(:)={{['Transport between ' fromComp ' and ' toComps{i}]}};
0140         else
0141             ssFiller(:)={{['Transport from ' fromComp ' to ' toComps{i}]}};
0142         end
0143         model.subSystems=[model.subSystems;ssFiller];
0144     end
0145     if isfield(model,'grRules')
0146         model.grRules=[model.grRules;filler];
0147     end
0148     if isfield(model,'rxnFrom')
0149         model.rxnFrom=[model.rxnFrom;filler];
0150     end
0151     if isfield(model,'rxnMiriams')
0152         model.rxnMiriams=[model.rxnMiriams;cell(nRxns,1)];
0153     end
0154     if isfield(model,'rxnComps')
0155         model.rxnComps=[model.rxnComps;ones(nRxns,1)];
0156         fprintf('NOTE: The added transport reactions will be added to the first compartment\n');
0157     end
0158     if isfield(model,'rxnGeneMat')
0159         model.rxnGeneMat=[model.rxnGeneMat;sparse(nRxns,numel(model.genes))];
0160     end
0161     if isfield(model,'rxnNotes')
0162         model.rxnNotes=[model.rxnNotes;filler];
0163     end
0164     if isfield(model,'rxnReferences')
0165         model.rxnReferences=[model.rxnReferences;filler];
0166     end
0167     if isfield(model,'rxnConfidenceScores')
0168         model.rxnConfidenceScores=[model.rxnConfidenceScores;ones(nRxns,1)];
0169     end
0170     if isfield(model,'rxnDeltaG')
0171         model.rxnDeltaG=[model.rxnDeltaG;zeros(nRxns,1)];
0172     end
0173     addedRxns = [addedRxns; addedRxnsID];
0174 end
0175 end