-
Notifications
You must be signed in to change notification settings - Fork 318
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'opencobra-develop' into develop
- Loading branch information
Showing
48 changed files
with
3,833 additions
and
1,496 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
795 changes: 357 additions & 438 deletions
795
...is/multiSpecies/microbiomeModelingToolbox/additionalAnalysis/calculateReactionAbundance.m
100755 → 100644
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
119 changes: 119 additions & 0 deletions
119
...ultiSpecies/microbiomeModelingToolbox/additionalAnalysis/fastCalculateReactionAbundance.m
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,119 @@ | ||
| function ReactionAbundance = fastCalculateReactionAbundance(abundancePath, modelPath, rxnsList, numWorkers) | ||
| % Part of the Microbiome Modeling Toolbox. This function calculates and | ||
| % plots the total abundance of reactions of interest in a given microbiome | ||
| % sample based on the strain-level composition. | ||
| % Reaction presence or absence in each strain is derived from the reaction content | ||
| % of the respective AGORA model. | ||
| % | ||
| % USAGE | ||
| % | ||
| % ReactionAbundance = fastCalculateReactionAbundance(abundancePath, modelPath, rxnsList, numWorkers) | ||
| % | ||
| % INPUTS: | ||
| % abundancePath: Path to the .csv file with the abundance data. | ||
| % Example: 'cobratoolbox/papers/018_microbiomeModelingToolbox/examples/normCoverage.csv' | ||
| % modelPath: Folder containing the strain-specific AGORA models | ||
| % OPTIONAL INPUTS: | ||
| % rxnsList: List of reactions for which the abundance | ||
| % should be calculated (if left empty: all | ||
| % reactions in all models) | ||
| % numWorkers: Number of workers used for parallel pool. If | ||
| % left empty, the parallel pool will not be | ||
| % started. Parallellization is recommended if | ||
| % all reactions are computed. | ||
| % | ||
| % OUTPUT: | ||
| % ReactionAbundance Table with total abundance for each microbiome | ||
| % and reaction | ||
| % | ||
| % .. Author: - Almut Heinken, 04/2021 | ||
|
|
||
| % read the csv file with the abundance data | ||
| abundance = readtable(abundancePath, 'ReadVariableNames', false); | ||
| abundance = table2cell(abundance); | ||
| if isnumeric(abundance{2, 1}) | ||
| abundance(:, 1) = []; | ||
| end | ||
|
|
||
| % load the models | ||
| for i = 2:size(abundance, 1) | ||
| model = readCbModel([modelPath filesep abundance{i, 1} '.mat']); | ||
| modelsList{i, 1} = model; | ||
| end | ||
|
|
||
| if ~exist('rxnsList', 'var') || isempty(rxnsList) % define reaction list if not entered | ||
| fprintf('No reaction list entered. Abundances will be calculated for all reactions in all models. \n') | ||
| % get model list from abundance input file | ||
| for i = 2:size(abundance, 1) | ||
| model = modelsList{i, 1}; | ||
| rxnsList = vertcat(model.rxns, rxnsList); | ||
| end | ||
| rxnsList = unique(rxnsList); | ||
| end | ||
|
|
||
| % load the models found in the individuals and extract which reactions are | ||
| % in which model | ||
| for i = 2:size(abundance, 1) | ||
| model = modelsList{i, 1}; | ||
| ReactionPresence{i, 1} = abundance{i, 1}; | ||
| for j = 1:length(rxnsList) | ||
| ReactionPresence{1, j + 1} = rxnsList{j}; | ||
| if ~isempty(find(ismember(model.rxns, rxnsList{j}))) | ||
| ReactionPresence{i, j + 1} = '1'; | ||
| else | ||
| ReactionPresence{i, j + 1} = '0'; | ||
| end | ||
| end | ||
| end | ||
| ReactionPresence{1,1}='Strains'; | ||
|
|
||
|
|
||
| % prepare table for the total abundance | ||
| ReactionAbundance = {}; | ||
| for i = 1:length(rxnsList) | ||
| ReactionAbundance{1, i + 1} = rxnsList{i}; | ||
| end | ||
| for i = 2:size(abundance, 2) | ||
| ReactionAbundance{i, 1} = abundance{1, i}; | ||
| end | ||
|
|
||
| % use parallel pool if workers specified as input | ||
| if exist('numWorkers', 'var') && numWorkers > 0 | ||
| poolobj = gcp('nocreate'); | ||
| if isempty(poolobj) | ||
| parpool(numWorkers) | ||
| end | ||
| end | ||
|
|
||
| clear abundance | ||
|
|
||
| totalAbun={}; | ||
| parfor i = 2:size(ReactionAbundance, 1) | ||
| i | ||
| % reload the file to avoid running out of memory | ||
| abundance = readtable(abundancePath, 'ReadVariableNames', false); | ||
| abundance = table2cell(abundance); | ||
| if isnumeric(abundance{2, 1}) | ||
| abundance(:, 1) = []; | ||
| end | ||
|
|
||
| % temporarily store reaction abundances | ||
| totalAbun{i} = zeros(length(rxnsList), 1); | ||
|
|
||
| for j = 2:size(abundance, 1) | ||
| % find all reactions present in the strain | ||
| presentRxns = find(strcmp(ReactionPresence(j,2:end),'1')); | ||
|
|
||
| for k = 1:length(presentRxns) | ||
| % summarize total abundance | ||
| totalAbun{i}(presentRxns(k),1) = totalAbun{i}(presentRxns(k),1) + str2double(abundance{j,i}); | ||
| end | ||
| end | ||
| end | ||
|
|
||
| % collect the temporarily stored abundances to put together the table | ||
| for i = 2:size(ReactionAbundance, 1) | ||
| ReactionAbundance(i,2:end) = num2cell(totalAbun{i}); | ||
| end | ||
|
|
||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
78 changes: 78 additions & 0 deletions
78
src/analysis/multiSpecies/microbiomeModelingToolbox/mgPipe/buildModelStorage.m
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,78 @@ | ||
| function [exch,modelStoragePath] = buildModelStorage(microbeNames,modPath) | ||
|
|
||
| currentDir=pwd; | ||
| mkdir('modelStorage') | ||
| cd('modelStorage') | ||
| modelStoragePath = pwd; | ||
|
|
||
| exch = {}; | ||
| for j = 1:size(microbeNames, 1) | ||
| model = readCbModel([modPath filesep microbeNames{j,1} '.mat']); | ||
| %exch = union(exch, model.mets(find(sum(model.S(:, strncmp('EX_', model.rxns, 3)), 2) ~= 0))); | ||
| exStruct = findSExRxnInd(model); | ||
| new_exch = findMetsFromRxns(model,model.rxns(exStruct.ExchRxnBool & ~exStruct.biomassBool)); | ||
| exch = union(exch,new_exch); | ||
| end | ||
|
|
||
| % get already built reconstructions | ||
| dInfo = dir(modelStoragePath); | ||
| modelList={dInfo.name}; | ||
| modelList=modelList'; | ||
| modelList=strrep(modelList,'.mat',''); | ||
| microbesNames=setdiff(microbeNames,modelList); | ||
|
|
||
|
|
||
| if length(microbesNames)>0 | ||
| %% create a new extracellular space [u] for microbes | ||
| for j = 1:size(microbeNames, 1) | ||
| model = readCbModel([modPath filesep microbeNames{j,1} '.mat']); | ||
| % temp fix | ||
| if isfield(model,'C') | ||
| model=rmfield(model,'C'); | ||
| model=rmfield(model,'d'); | ||
| end | ||
| % | ||
|
|
||
| % removing possible constraints of the bacs | ||
| selExc = findExcRxns(model); | ||
| Reactions2 = model.rxns(find(selExc)); | ||
| allex = Reactions2(strmatch('EX', Reactions2)); | ||
| biomass = allex(find(strncmp(allex,'bio',3))); | ||
| finrex = setdiff(allex, biomass); | ||
| model = changeRxnBounds(model, finrex, -1000, 'l'); | ||
| model = changeRxnBounds(model, finrex, 1000, 'u'); | ||
|
|
||
| % removing blocked reactions from the bacs | ||
| %BlockedRxns = identifyFastBlockedRxns(model,model.rxns, printLevel); | ||
| %model= removeRxns(model, BlockedRxns); | ||
| %BlockedReaction = findBlockedReaction(model,'L2') | ||
|
|
||
| model = convertOldStyleModel(model); | ||
| exmod = model.rxns(strncmp('EX_', model.rxns, 3)); % find exchange reactions | ||
| eMets = model.mets(~cellfun(@isempty, strfind(model.mets, '[e]'))); % exchanged metabolites | ||
| dummyMicEU = createModel(); | ||
| %dummyMicEU = makeDummyModel(2 * size(eMets, 1), size(eMets, 1)); | ||
| dummyMicEUmets = [strcat(strcat(microbeNames{j, 1}, '_'), regexprep(eMets, '\[e\]', '\[u\]')); regexprep(eMets, '\[e\]', '\[u\]')]; | ||
| dummyMicEU = addMultipleMetabolites(dummyMicEU,dummyMicEUmets); | ||
| nMets = numel(eMets); | ||
| S = [speye(nMets);-speye(nMets)]; | ||
| lbs = repmat(-1000,nMets,1); | ||
| ubs = repmat(1000,nMets,1); | ||
| names = strcat(strcat(microbeNames{j, 1}, '_'), 'IEX_', regexprep(eMets, '\[e\]', '\[u\]'), 'tr'); | ||
| dummyMicEU = addMultipleReactions(dummyMicEU,names,dummyMicEUmets,S,'lb',lbs,'ub',ubs); | ||
| model = removeRxns(model, exmod); | ||
| model.rxns = strcat(strcat(microbeNames{j, 1}, '_'), model.rxns); | ||
| model.mets = strcat(strcat(microbeNames{j, 1}, '_'), regexprep(model.mets, '\[e\]', '\[u\]')); % replace [e] with [u] | ||
| [model] = mergeTwoModels(dummyMicEU, model, 2, false, false); | ||
|
|
||
| %finish up by A: removing duplicate reactions | ||
| %We will lose information here, but we will just remove the duplicates. | ||
| [model,rxnToRemove,rxnToKeep]= checkDuplicateRxn(model,'S',1,0,1); | ||
|
|
||
| writeCbModel(model,'format','mat','fileName',[microbeNames{j,1} '.mat']); % store model | ||
| end | ||
| end | ||
|
|
||
| cd(currentDir) | ||
|
|
||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.