Merge pull request #2371 from opencobra/develop

Develop

src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/analyzeMgPipeResults.m

@@ -112,9 +112,9 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
             % Print the results as a text file
             filename = strrep(fileList{i},'.csv','');
             filename = strrep(filename,'.txt','');
-            writetable(cell2table(Statistics),[statPath filesep filename '_' sampleGroupHeaders{j} '_Statistics'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+            writetable(cell2table(Statistics),[statPath filesep filename '_' sampleGroupHeaders{j} '_Statistics.csv'],'WriteVariableNames',false);
             if size(significantFeatures,2)>1
-                writetable(cell2table(significantFeatures),[statPath filesep filename '_' sampleGroupHeaders{j} '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+                writetable(cell2table(significantFeatures),[statPath filesep filename '_' sampleGroupHeaders{j} '_SignificantFeatures.csv'],'WriteVariableNames',false);
             end
 
             % create plots
@@ -133,9 +133,9 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
         % Print the results as a text file
         filename = strrep(fileList{i},'.csv','');
         filename = strrep(filename,'.txt','');
-        writetable(cell2table(Statistics),[statPath filesep filename '_Statistics'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+        writetable(cell2table(Statistics),[statPath filesep filename '_Statistics.csv'],'WriteVariableNames',false);
         if size(significantFeatures,2)>1
-            writetable(cell2table(significantFeatures),[statPath filesep filename '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
+            writetable(cell2table(significantFeatures),[statPath filesep filename '_SignificantFeatures.csv'],'WriteVariableNames',false);
         end
 
         % create violin plots

src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/makeViolinPlots.m

@@ -54,9 +54,16 @@ function makeViolinPlots(sampleData, sampleInformation, varargin)
 cnt=1;
 delArray=[];
 for i=2:size(sampleData,1)
-    if sum(str2double(sampleData(i,2:end)))<0.0001
-        delArray(cnt,1)=i;
-        cnt=cnt+1;
+    if contains(version,'R202') % for MATLAB 2020a or newer
+        if sum(cell2mat(sampleData(i,2:end)))<0.0001
+            delArray(cnt,1)=i;
+            cnt=cnt+1;
+        end
+    else
+        if sum(str2double(sampleData(i,2:end)))<0.0001
+            delArray(cnt,1)=i;
+            cnt=cnt+1;
+        end
     end
 end
 sampleData(delArray,:)=[];

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

@@ -50,6 +50,8 @@
 % define human-derived metabolites present in the gut: primary bile acids, amines, mucins, host glycans
 if includeHumanMets
     HumanMets={'gchola','-10';'tdchola','-10';'tchola','-10';'dgchol','-10';'34dhphe','-10';'5htrp','-10';'Lkynr','-10';'f1a','-1';'gncore1','-1';'gncore2','-1';'dsT_antigen','-1';'sTn_antigen','-1';'core8','-1';'core7','-1';'core5','-1';'core4','-1';'ha','-1';'cspg_a','-1';'cspg_b','-1';'cspg_c','-1';'cspg_d','-1';'cspg_e','-1';'hspg','-1'};
+else
+    HumanMets={};
 end
 
 % get all exchanges that can carry flux in at least one model on the given
@@ -94,11 +96,14 @@
             end
         end
         model = useDiet(model, diet,0);
-
+        
         if includeHumanMets
-            % add the human metabolites
+            % add the human metabolites if not already included
+            % in the diet
             for l=1:length(HumanMets)
-                model=changeRxnBounds(model,strcat('EX_',HumanMets{l},'(e)'),str2num(HumanMets{l,2}),'l');
+                if isempty(find(strcmp(diet(:,1),['EX_',HumanMets{l},'(e)'])))
+                    model=changeRxnBounds(model,['EX_',HumanMets{l},'(e)'],str2num(HumanMets{l,2}),'l');
+                end
             end
         end
     end
@@ -174,11 +179,14 @@
                 end
             end
             model = useDiet(model, diet,0);
-
+            
             if includeHumanMets
-                % add the human metabolites
+                % add the human metabolites if not already included
+                % in the diet
                 for l=1:length(HumanMets)
-                    model=changeRxnBounds(model,strcat('EX_',HumanMets{l},'(e)'),str2num(HumanMets{l,2}),'l');
+                    if isempty(find(strcmp(diet(:,1),['EX_',HumanMets{l},'(e)'])))
+                        model=changeRxnBounds(model,['EX_',HumanMets{l},'(e)'],str2num(HumanMets{l,2}),'l');
+                    end
                 end
             end
         end

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

@@ -1,4 +1,4 @@
-function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorkers)
+function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorkers, builtTaxa)
 % This function creates pan-models for all unique taxa (e.g., species)
 % included in the AGORA resource. If reconstructions of multiple strains
 % in a given taxon are present, the reactions in these reconstructions will
@@ -23,16 +23,22 @@ function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorker
 %                  created pan-models will be stored in. Must end with a file separator.
 %    taxonLevel    String with desired taxonomical level of the pan-models.
 %                  Allowed inputs are 'Species','Genus','Family','Order', 'Class','Phylum'.
-%    agoraVersion  Version of AGORA that will be used (allowed inputs: 'AGORA', 'AGORA2')
+%    agoraVersion  Version of AGORA that will be used (allowed inputs: 'AGORA', 'AGORA2', 
+%                  alternatively: path to custom table with reconstruction information)
 %
 % OPTIONAL INPUTS
 %    numWorkers    Number of workers for parallel pool (default: no pool)
+%    builtTaxa     Names of taxa in table that will be built (default:
+%                  all). Need to be entered as a cell array of strings with names written
+%                  exactly as in the corresponding column in the table.
 %
 % .. Authors
 %       - Stefania Magnusdottir, 2016
 %       - Almut Heinken, 06/2018: adapted to function.
 %       - Almut Heinken, 03/2021: enabled parallelization
 %       - Almut Heinken, 05/2024: enforced specifying AGORA version
+%       - Almut Heinken, 10/2024: allowed specifying input table and taxa
+%                                 to create
 
 tol = 1e-5;
 
@@ -58,13 +64,13 @@ function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorker
 solver = CBT_LP_SOLVER;
 environment = getEnvironment();
 
-% load the file with information on AGORA/AGORA2
+% load the file with information on AGORA/AGORA2 or custom table
 if strcmp(agoraVersion,'AGORA')
     infoFile = readInputTableForPipeline('AGORA_infoFile.xlsx');
 elseif strcmp(agoraVersion,'AGORA2')
     infoFile = readInputTableForPipeline('AGORA2_infoFile.xlsx');
 else
-    error('Incorrect input for agoraVersion! Allowed inputs: AGORA, AGORA2')
+    infoFile = readInputTableForPipeline(agoraVersion);
 end
 
 % get the reaction and metabolite database
@@ -78,6 +84,12 @@ function createPanModels(agoraPath, panPath, taxonLevel, agoraVersion, numWorker
 allTaxa(strncmp(allTaxa, 'unclassified',  12)) = [];
 allTaxa(strcmp(allTaxa, '')) = [];
 
+% reduce to specified list of taxa to build if entered
+if nargin >5
+    [~,I] = setdiff(allTaxa,builtTaxa);
+    allTaxa(I,:) = [];
+end
+
 % Remove models that have already been assembled from the list of models to create
 dInfo = dir(panPath);
 dInfo = dInfo(~[dInfo.isdir]);

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

@@ -79,7 +79,11 @@
 allDietExch = regexprep(allDietExch,'\[fe\]','\[d\]');
 
 % define human-derived metabolites present in the gut: primary bile acids, amines, mucins, host glycans
-HumanMets={'gchola','-10';'tdchola','-10';'tchola','-10';'dgchol','-10';'34dhphe','-10';'5htrp','-10';'Lkynr','-10';'f1a','-1';'gncore1','-1';'gncore2','-1';'dsT_antigen','-1';'sTn_antigen','-1';'core8','-1';'core7','-1';'core5','-1';'core4','-1';'ha','-1';'cspg_a','-1';'cspg_b','-1';'cspg_c','-1';'cspg_d','-1';'cspg_e','-1';'hspg','-1'};
+if includeHumanMets
+    HumanMets={'gchola','-10';'tdchola','-10';'tchola','-10';'dgchol','-10';'34dhphe','-10';'5htrp','-10';'Lkynr','-10';'f1a','-1';'gncore1','-1';'gncore2','-1';'dsT_antigen','-1';'sTn_antigen','-1';'core8','-1';'core7','-1';'core5','-1';'core4','-1';'ha','-1';'cspg_a','-1';'cspg_b','-1';'cspg_c','-1';'cspg_d','-1';'cspg_e','-1';'hspg','-1'};
+else
+    HumanMets = {};
+end
 
 %% start the simulations
 
@@ -267,7 +271,7 @@
                     model=changeRxnBounds(model,['Host_' hostBiomassRxn],0.001,'l');
                     model=changeRxnBounds(model,['Host_' hostBiomassRxn],hostBiomassRxnFlux,'u');
                 end
-
+                
                 solution_allOpen = optimizeCbModel(model);
                 % solution_allOpen=solveCobraLPCPLEX(model,2,0,0,[],0);
                 if solution_allOpen.stat==0
@@ -281,13 +285,13 @@
                     FecalRxn = AllRxn(FecalInd);
                     FecalRxn=setdiff(FecalRxn,'EX_microbeBiomass[fe]','stable');
                     DietRxn = AllRxn(DietInd);
-
+                    
                     %% computing fluxes on the rich diet
                     if rDiet==1 && computeProfiles
                         % remove exchanges that cannot carry flux
                         FecalRxn=intersect(FecalRxn,allFecalExch);
                         DietRxn=intersect(DietRxn,allDietExch);
-
+                        
                         [minFlux,maxFlux]=guidedSim(model,FecalRxn);
                         minFluxFecal = minFlux;
                         maxFluxFecal = maxFlux;
@@ -304,32 +308,35 @@
                             netUptakeTmp{k}{1}{index,3} = minFluxFecal(i,1);
                         end
                     end
-
+                    
                     %% Computing fluxes on the input diet
-
+                    
                     % remove exchanges that cannot carry flux
                     FecalRxn=intersect(FecalRxn,allFecalExch);
                     DietRxn=intersect(DietRxn,allDietExch);
-
+                    
                     model_sd=model;
                     if adaptMedium
                         [diet] = adaptVMHDietToAGORA(loadDiet,'Microbiota');
                     else
                         diet = readInputTableForPipeline(loadDiet);  % load the text file with the diet
-
-                        for j = 2:length(diet)
+                        
+                        for j = 1:length(diet)
                             diet{j, 2} = num2str(-(diet{j, 2}));
                         end
                     end
                     [model_sd] = useDiet(model_sd, diet,0);
-
+                    
                     if includeHumanMets
-                        % add the human metabolites
+                        % add the human metabolites if not already included
+                        % in the diet
                         for l=1:length(HumanMets)
-                            model_sd=changeRxnBounds(model_sd,strcat('Diet_EX_',HumanMets{l},'[d]'),str2num(HumanMets{l,2}),'l');
+                            if isempty(find(strcmp(diet(:,1),['Diet_EX_',HumanMets{l},'[d]'])))
+                                model_sd=changeRxnBounds(model_sd,['Diet_EX_',HumanMets{l},'[d]'],str2num(HumanMets{l,2}),'l');
+                            end
                         end
                     end
-
+                    
                     solution_sDiet=optimizeCbModel(model_sd);
                     % solution_sDiet=solveCobraLPCPLEX(model_sd,2,0,0,[],0);
                     growthRatesTmp{k}{2}=solution_sDiet.f;
@@ -356,64 +363,64 @@
                                 netUptakeTmp{k}{2}{index,3} = minFluxFecal(i,1);
                             end
                         end
-
-                        microbiota_model=model_sd;
-                        parsave([resPath filesep 'Diet' filesep 'microbiota_model_diet_' sampleID '.mat'],microbiota_model)
-
-                        %% Using personalized diet not documented in MgPipe and bug checked yet!!!!
-
-                        if pDiet==1
-                            model_pd=model;
-                            [Numbers, Strings] = xlsread(strcat(abundancepath,fileNameDiets));
-                            % diet exchange reactions
-                            DietNames = Strings(2:end,1);
-                            % Diet exchanges for all individuals
-                            Diets(:,k) = cellstr(num2str((Numbers(1:end,k))));
-                            Dietexchanges = {DietNames{:,1} ; Diets{:,k}}';
-                            Dietexchanges = regexprep(Dietexchanges,'EX_','Diet_EX_');
-                            Dietexchanges = regexprep(Dietexchanges,'\(e\)','\[d\]');
-
-                            model_pd = setDietConstraints(model_pd,Dietexchanges);
-
-                            if includeHumanMets
-                                % add the human metabolites
-                                for l=1:length(HumanMets)
-                                    model_pd=changeRxnBounds(model_pd,strcat('Diet_EX_',HumanMets{l},'[d]'),str2num(HumanMets{l,2}),'l');
-                                end
+                    end
+
+                    microbiota_model=model_sd;
+                    parsave([resPath filesep 'Diet' filesep 'microbiota_model_diet_' sampleID '.mat'],microbiota_model)
+
+                    %% Using personalized diet not documented in MgPipe and bug checked yet!!!!
+
+                    if pDiet==1
+                        model_pd=model;
+                        [Numbers, Strings] = xlsread(strcat(abundancepath,fileNameDiets));
+                        % diet exchange reactions
+                        DietNames = Strings(2:end,1);
+                        % Diet exchanges for all individuals
+                        Diets(:,k) = cellstr(num2str((Numbers(1:end,k))));
+                        Dietexchanges = {DietNames{:,1} ; Diets{:,k}}';
+                        Dietexchanges = regexprep(Dietexchanges,'EX_','Diet_EX_');
+                        Dietexchanges = regexprep(Dietexchanges,'\(e\)','\[d\]');
+
+                        model_pd = setDietConstraints(model_pd,Dietexchanges);
+
+                        if includeHumanMets
+                            % add the human metabolites
+                            for l=1:length(HumanMets)
+                                model_pd=changeRxnBounds(model_pd,strcat('Diet_EX_',HumanMets{l},'[d]'),str2num(HumanMets{l,2}),'l');
                             end
-
-                            solution_pdiet=optimizeCbModel(model_pd);
-                            %solution_pdiet=solveCobraLPCPLEX(model_pd,2,0,0,[],0);
-                            growthRatesTmp{k}{3}=solution_pdiet.f;
-                            if solution_pdiet.stat==0
-                                warning('growthRates detected one or more infeasible models. Please check infeasModels object !')
-                                infeasModelsTmp{k} = model.name;
-                                netProductionTmp{k}{3} = {};
-                                netUptakeTmp{k}{3} = {};
-                            else
-                                if computeProfiles
-                                    [minFlux,maxFlux]=guidedSim(model_pd,FecalRxn);
-                                    minFluxFecal = minFlux;
-                                    maxFluxFecal = maxFlux;
-                                    [minFlux,maxFlux]=guidedSim(model_pd,DietRxn);
-                                    minFluxDiet = minFlux;
-                                    maxFluxDiet = maxFlux;
-                                    netProductionTmp{k}{3}=exchanges;
-                                    netUptakeTmp{k}{3}=exchanges;
-                                    for i =1:length(FecalRxn)
-                                        [truefalse, index] = ismember(FecalRxn(i), exchanges);
-                                        netProductionTmp{k}{3}{index,2} = minFluxDiet(i,1);
-                                        netProductionTmp{k}{3}{index,3} = maxFluxFecal(i,1);
-                                        netUptakeTmp{k}{3}{index,2} = maxFluxDiet(i,1);
-                                        netUptakeTmp{k}{3}{index,3} = minFluxFecal(i,1);
-                                    end
+                        end
+
+                        solution_pdiet=optimizeCbModel(model_pd);
+                        %solution_pdiet=solveCobraLPCPLEX(model_pd,2,0,0,[],0);
+                        growthRatesTmp{k}{3}=solution_pdiet.f;
+                        if solution_pdiet.stat==0
+                            warning('growthRates detected one or more infeasible models. Please check infeasModels object !')
+                            infeasModelsTmp{k} = model.name;
+                            netProductionTmp{k}{3} = {};
+                            netUptakeTmp{k}{3} = {};
+                        else
+                            if computeProfiles
+                                [minFlux,maxFlux]=guidedSim(model_pd,FecalRxn);
+                                minFluxFecal = minFlux;
+                                maxFluxFecal = maxFlux;
+                                [minFlux,maxFlux]=guidedSim(model_pd,DietRxn);
+                                minFluxDiet = minFlux;
+                                maxFluxDiet = maxFlux;
+                                netProductionTmp{k}{3}=exchanges;
+                                netUptakeTmp{k}{3}=exchanges;
+                                for i =1:length(FecalRxn)
+                                    [truefalse, index] = ismember(FecalRxn(i), exchanges);
+                                    netProductionTmp{k}{3}{index,2} = minFluxDiet(i,1);
+                                    netProductionTmp{k}{3}{index,3} = maxFluxFecal(i,1);
+                                    netUptakeTmp{k}{3}{index,2} = maxFluxDiet(i,1);
+                                    netUptakeTmp{k}{3}{index,3} = minFluxFecal(i,1);
                                 end
-
-                                % save the model with personalized diet
-                                microbiota_model=model_pd;
-                                mkdir(strcat(resPath,'Personalized'))
-                                parsave([resPath filesep 'Personalized' filesep 'microbiota_model_pDiet_' sampleID '.mat'],microbiota_model)
                             end
+
+                            % save the model with personalized diet
+                            microbiota_model=model_pd;
+                            mkdir(strcat(resPath,'Personalized'))
+                            parsave([resPath filesep 'Personalized' filesep 'microbiota_model_pDiet_' sampleID '.mat'],microbiota_model)
                         end
                     end
                 end
@@ -435,14 +442,14 @@
                 end
             end
             if ~isempty(growthRatesTmp{k})
-                 if ~isempty(growthRatesTmp{k}{1})
+                if ~isempty(growthRatesTmp{k}{1})
                     growthRates{k+1,2} = growthRatesTmp{k}{1};
                     growthRates{k+1,3} = growthRatesTmp{k}{2};
                     if length(growthRatesTmp{k})>2
                         growthRates{1,4} = 'Personalized diet';
                         growthRates{k+1,4} = growthRatesTmp{k}{3};
                     end
-                 end
+                end
             end
             if ~isempty(infeasModelsTmp) && k <= length(infeasModelsTmp)
                 infeasModels{k,1} = infeasModelsTmp{k};