Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixed some bugs in DEMETER/mgPipe #2350

Merged
merged 6 commits into from
Nov 7, 2024
Merged

Fixed some bugs in DEMETER/mgPipe #2350

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
cd75070
Fixed some bugs in DEMETER/mgPipe
almut-heinken Oct 8, 2024
b3492e4
Fixed some bugs in DEMETER/mgPipe
almut-heinken Oct 10, 2024
b9daa66
Fixed some bugs in DEMETER/mgPipe
almut-heinken Oct 11, 2024
9cf77b5
Fixed some bugs in DEMETER/mgPipe
almut-heinken Oct 14, 2024
a369b99
Fixed some bugs in DEMETER/mgPipe
almut-heinken Oct 16, 2024
15b8766
added distributed models for testCreatePanModels
almut-heinken Nov 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions ...analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/analyzeMgPipeResults.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down
13 changes: 10 additions & 3 deletions src/analysis/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/makeViolinPlots.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,:)=[];
Expand Down
20 changes: 14 additions & 6 deletions src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/buildModelStorage.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
20 changes: 16 additions & 4 deletions src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/createPanModels.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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;

Expand All @@ -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
Expand All @@ -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]);
Expand Down
147 changes: 77 additions & 70 deletions src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/microbiotaModelSimulator.m
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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
Expand All @@ -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};
Expand Down
Loading
Loading