Merge pull request #1644 from opencobra/develop

Merger of PSCM toolbox described in PMID 32463598

external/analysis/StoichTools/molweight.m

@@ -45,6 +45,8 @@
     if isempty(ppds)
         element = @(str, atomicnumber, amu) amu;
 
+        ppds.A  = element('R group',      0,   NaN);%Added by Ronan Fleming
+
         ppds.M  = element('any metal',      0,   NaN);
         ppds.X  = element('any halogen',    0,   NaN);
 

external/analysis/StoichTools/parse_formula.m

@@ -65,6 +65,10 @@
 %   6. The bare electron e- is used in balancing chemical half reactions.
 %
 %   7. Error messages are generated for invalid fomulas
+%      Change for use with the cobra toolbox
+%      A is accepted as the symbol for an element that could
+%      represent an R group in an atom mapping scenario for an incompletely
+%      defined metabolite.
 %
 %   8. str can be a cell array of chemical formula. The results is a
 %      structure array. The elements of the output structure array are in
@@ -162,7 +166,9 @@
     persistent grex;  % Regexp pattern to match groups
 
     if isempty(srex) || isempty(grex)
-        srex = ['(A[lrsgutcm]|B[eraik]?|C[laroudsemf]?|D[y]?|E[urs]|', ...
+        %this was the original line
+        %srex = ['(A[lrsgutcm]|B[eraik]?|C[laroudsemf]?|D[y]?|E[urs]|',...
+        srex = ['(A[lrsgutcm]?|B[eraik]?|C[laroudsemf]?|D[y]?|E[urs]|', ...
                 'F[erm]?|G[aed]|H[eofgas]?|I[nr]?|Kr?|L[iaur]|', ...
                 'M[gnodt]?|N[eaibdpos]?|Os?|P[drmtboau]?|R[buhenaf]|', ...
                 'S[icernbmg]?|T[icebmalh]?|U|V|W|X[e]?|Yb?|Z[nr])', ...

src/analysis/FBA/optimizeCbModel.m

@@ -432,6 +432,7 @@
         solution.full   = solutionL0.x;
         solution.dual   = [];
         solution.rcost  = [];
+        solution.slack  = [];
 
     elseif length(minNorm)> 1 || minNorm > 0
         %THIS SECTION BELOW ASSUMES WRONGLY THAT c HAVE ONLY ONE NONZERO SO I

src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/addMicrobeCommunityBiomass.m

@@ -1,44 +1,50 @@
-function model = addMicrobeCommunityBiomass(model, microbeNames, abundances)
-% Adds a community biomass reaction to a model structure with multiple
-% microbes based on their relative abundances. If no abundance values are
-% provided, all n microbes get equal weights (1/n). Assumes a lumen
-% compartment [u] and a fecal secretion comparment [fe]. Creates a community
-% biomass metabolite 'microbeBiomass' that is secreted from [u] to [fe] and
-% exchanged from fecal compartment.
-%
-% USAGE:
-%
-%   model = addMicrobeCommunityBiomass(model, microbeNames, abundances)
-%
-% INPUTS:
-%   model:           COBRA model structure with n joined microbes with biomass
-%                    metabolites 'Microbe_biomass[c]'.
-%   microbeNames:    nx1 cell array of n unique strings that represent
-%                    each microbe in the model.
-%
-% OPTIONAL INPUT:
-%   abundances:      nx1 vector with the relative abundance of each microbe.
-%
-% OUTPUT:
-%   model:           COBRA model structure
-%
-% .. Author: Stefania Magnusdottir June 2016
-
-dummy = createModel(); %makeDummyModel(length(microbeNames) + 2, 3);
-
-mets = [strcat(microbeNames, '_biomass[c]'); 'microbeBiomass[u]'; 'microbeBiomass[fe]'];
-rxns = {'communityBiomass'; 'UFEt_microbeBiomass'; 'EX_microbeBiomass[fe]'};
-if ~exist('abundances','var') || isempty(abundances)    
-    S = [-ones(size(microbeNames)) / length(microbeNames); 1; 0];    
-else
-    S = [-abundances; 1; 0];
-end
-S(end-1:end,2) = [-1; 1];
-S(end,3) = -1;
-% three reactions are added
-lb = [0;0;-1000];
-ub = ones(3,1) * 1000;
-dummy = addMultipleMetabolites(dummy,mets);
-dummy = addMultipleReactions(dummy,rxns,mets,S,'lb',lb,'ub',ub);
-% join models
-model = mergeTwoModels(dummy, model, 2, 0);
+function model = addMicrobeCommunityBiomass(model, microbeNames, abundances)
+% Adds a community biomass reaction to a model structure with multiple
+% microbes based on their relative abundances. If no abundance values are
+% provided, all n microbes get equal weights (1/n). Assumes a lumen
+% compartment [u] and a fecal secretion comparment [fe]. Creates a community
+% biomass metabolite 'microbeBiomass' that is secreted from [u] to [fe] and
+% exchanged from fecal compartment.
+%
+% USAGE:
+%
+%   model = addMicrobeCommunityBiomass(model, microbeNames, abundances)
+%
+% INPUTS:
+%   model:           COBRA model structure with n joined microbes with biomass
+%                    metabolites 'Microbe_biomass[c]'.
+%   microbeNames:    nx1 cell array of n unique strings that represent
+%                    each microbe in the model.
+%
+% OPTIONAL INPUT:
+%   abundances:      nx1 vector with the relative abundance of each microbe.
+%
+% OUTPUT:
+%   model:           COBRA model structure
+%
+% .. Author: Stefania Magnusdottir June 2016
+
+dummy = createModel(); %makeDummyModel(length(microbeNames) + 2, 3);
+
+mets = [strcat(microbeNames, '_biomass[c]'); 'microbeBiomass[u]'; 'microbeBiomass[fe]'];
+rxns = {'communityBiomass'; 'UFEt_microbeBiomass'; 'EX_microbeBiomass[fe]'};
+if ~exist('abundances','var') || isempty(abundances)
+    S = [-ones(size(microbeNames)) / length(microbeNames); 1; 0];
+else
+    S = [-abundances; 1; 0];
+end
+S(end-1:end,2) = [-1; 1];
+S(end,3) = -1;
+% three reactions are added
+lb = [0;0;-1000];
+ub = ones(3,1) * 1000;
+dummy = addMultipleMetabolites(dummy,mets);
+dummy = addMultipleReactions(dummy,rxns,mets,S,'lb',lb,'ub',ub);
+% join models
+%model = mergeTwoModels(dummy, model, 2, 0);
+% modification made by IT, May 2020
+% alternative and faster way to add the 3 reactions
+a = printRxnFormula(dummy);
+for i = 1 : length(a)
+    model = addReaction(model, dummy.rxns{i},a{i});
+end

src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/initMgPipe.m

@@ -39,9 +39,7 @@
 %               - Almut Heinken 02/2020: removed unnecessary outputs
 %               - Almut Heinken 08/2020: added extra inputs and changed to
 %                                        varargin input
-
-global CBTDIR
-
+%
 % Define default input parameters if not specified
 parser = inputParser();
 parser.addRequired('modPath', @ischar);
@@ -86,6 +84,13 @@
 lowerBMBound = parser.Results.lowerBMBound;
 repeatSim = parser.Results.repeatSim;
 
+global CBT_LP_SOLVER
+if isempty(CBT_LP_SOLVER)
+    initCobraToolbox
+end
+
+global CBTDIR
+
 % set optional variables
 if ~exist('resPath', 'var') || ~exist(resPath, 'dir')
     resPath = [CBTDIR filesep '.tmp'];

src/analysis/thermo/molFiles/kegg2mol.m

@@ -147,8 +147,14 @@ function kegg2mol(cid, molfileDir, mets, takeMajorMS, pH, takeMajorTaut)
             majorTautOption = 'false';
         end
 
-        status = system(['cxcalc -o ' molfileDir met '.mol majorms -H ' num2str(pH) ' -f mol -M ' majorTautOption ' tmp.mol']); % Call ChemAxon's calculator plugin (cxcalc) to compute major microspecies
-
+        if ismac
+            %expecting the default installation location
+           status = system(['/Applications/MarvinSuite/bin/cxcalc -o ' molfileDir met '.mol majorms -H ' num2str(pH) ' -f mol -M ' majorTautOption ' tmp.mol']); % Call ChemAxon's calculator plugin (cxcalc) to compute major microspecies
+        else
+            %user must make sure cxcalc is on their $PATH
+           status = system(['cxcalc -o ' molfileDir met '.mol majorms -H ' num2str(pH) ' -f mol -M ' majorTautOption ' tmp.mol']); % Call ChemAxon's calculator plugin (cxcalc) to compute major microspecies
+        end
+
         if status ~= 0
             nomol = [nomol; {met}];
         end

src/analysis/topology/conservedMoieties/generateMoietyMoleculeNetwork.m

@@ -0,0 +1,74 @@
+function MMN = generateMoietyMoleculeNetwork(model,L,moietyFormulae,mbool)
+%Generate a structure to represent a network that associates moieties to molecules
+
+[nMet,~]=size(model.S);
+[nMoieties, nMappedMet] = size(L);
+
+if ~exist('mbool','var')
+    if nMet==nMappedMet
+        mbool=true(nMet,1);
+    else
+        if nMet==nMoieties
+            L=L';
+            [nMoieties, nMappedMet] = size(L);
+        else
+            error('mbool must be provided if the number of columns of L is different from the number of rows of model.S')
+        end
+    end
+else
+    if nMappedMet~=nnz(mbool)
+        if nMoieties==nnz(mbool)
+            L=L';
+            [nMoieties, nMappedMet] = size(L);
+        else
+            error('mbool must have the same number of nonzeros as the number of columns of L')
+        end
+    end
+end
+
+MMN=struct();
+
+MMN.L = L;
+
+%order of the numbers matches the order of the rows of L
+moietyID = zeros(nMoieties,1);
+for i=1:nMoieties
+    moietyID(i) = i;
+end
+
+moietyFormulae = hillformula(moietyFormulae);
+
+isotopeAbundance = 0; %use polyisotopic inexact mass i.e. uses all isotopes of each element weighted by natural abundance 
+generalFormula = 1; %NaN for unknown elements
+[moietyMasses, knownMasses, unknownElements, Ematrix, elements] = getMolecularMass(moietyFormulae,isotopeAbundance,generalFormula);
+
+NumMolecules = sum(L,2);
+
+%table of moiety properties
+MMN.moi = table(moietyID,moietyFormulae,moietyMasses,NumMolecules,'VariableNames',{'Name','Formula','Mass','NumMolecules'});
+
+if 0
+    %     Error using parse_formula>parse_formula_ (line 197)
+    %     Could not parse formula:
+    %     C30H56NO12FULLRCO
+    %     ^^^
+    molFormulae = hillformula(model.metFormulas(mbool));
+else
+    molFormulae = model.metFormulas(mbool);
+end
+
+%order of the numbers matches the order of the columns of L
+molID = zeros(nMappedMet,1);
+for i=1:nMappedMet
+    molID(i) = i;
+end
+
+[molecularMasses, knownMasses, unknownElements, Ematrix, elements] = getMolecularMass(model.metFormulas(mbool),isotopeAbundance,generalFormula);
+
+NumMoieties = sum(L,1)';
+
+%table of molecule properties
+MMN.mol = table(molID,model.mets(mbool),molFormulae,molecularMasses,NumMoieties,'VariableNames',{'Name','Mets','Formula','Mass','NumMoieties'});
+
+end
+