Develop

.gitmodules

@@ -75,3 +75,6 @@
 [submodule "external/base/io/xlread"]
 	path = external/base/io/xlread
 	url = https://github.com/tpfau/xlread
+[submodule "external/visualization/MatGPT"]
+	path = external/visualization/MatGPT
+	url = [email protected]:toshiakit/MatGPT.git

src/analysis/FVA/fastFVA/mex/generateMexFastFVA.m

@@ -117,7 +117,7 @@ function generateMexFastFVA(rootPathCPLEX, printLevel)
     eval(['mex -setup c' compVerboseMode]);
 
     % define the name of the source code
-    filename = [CBTDIR filesep 'external' filesep 'analysis' filesep 'fastFVAmex' filesep 'cplexFVA.c'];
+    filename = [CBTDIR filesep 'src' filesep 'analysis' filesep 'FVA' filesep 'fastFVA' filesep 'mex' filesep 'cplexFVA.c'];
 
     % generation of MEX string with compiler options
     if isunix

src/analysis/calculateFluxShifts/calculateFluxShifts.m

@@ -1,6 +1,6 @@
 function calculateFluxShifts(source, target)
 % Compares reaction flux values (min, max) between two models or more models
-% in resultsPostOptimization\contextSpecificModels folder and writes the
+% in resultsPostOptimization/contextSpecificModels folder and writes the
 % comparison results ('up', 'down', 1) in a new column in the target model
 %
 % USAGE:
@@ -16,22 +16,23 @@ function calculateFluxShifts(source, target)
 % .. Authors:
 %       - Kristina Grausa 05/16/2022
 %       - Kristina Grausa 08/23/2022 - standard header and formatting
+%       - Farid Zare      20/11/2023 - Repository addresses are corrected, standard variable names
 
     % Get target model/-s
     destTarget = {};
     if string(target) == "All"
-        destTarget = string(strcat('resultsPostOptimization\contextSpecificModels\*.xls'));
+        destTarget = string(strcat('resultsPostOptimization/contextSpecificModels/*.xls'));
     else
-        destTarget = string(strcat('resultsPostOptimization\contextSpecificModels\*',char(target),'*.xls'));
+        destTarget = string(strcat('resultsPostOptimization/contextSpecificModels/*',char(target),'*.xls'));
     end
-    T = dir(destTarget);
+    xlsFilesInfo = dir(destTarget);
 
     % Source model
-    destSource = string(strcat('resultsPostOptimization\contextSpecificModels\*',char(source),'*.xls'));
+    destSource = string(strcat('resultsPostOptimization/contextSpecificModels/*',char(source),'*.xls'));
     S = dir(destSource);
     sheetName = 'Reaction List';
-    filename = strcat(S.folder, '\', S.name);
-    dataSource=readtable(filename,'Sheet',sheetName);
+    filename = strcat(S.folder, '/', S.name);
+    dataSource = readtable(filename,'Sheet',sheetName);
 
     % Get minFlux and MaxFlux columns from source model
     minFluxSource = dataSource.MinFlux;
@@ -47,59 +48,59 @@ function calculateFluxShifts(source, target)
     end
 
     % Get min and max flux data from target model file
-    for i=1:1:length(T)
-        if ~contains(T(i).name,char(source))
+    for i=1:length(xlsFilesInfo)
 
-            filename = strcat(T(i).folder, '\', T(i).name);
+        if ~contains(xlsFilesInfo(i).name,char(source))
+            filename = strcat(xlsFilesInfo(i).folder, '/', xlsFilesInfo(i).name);
             data=readtable(filename,'Sheet',sheetName);
 
             try
             % Get minFlux and MaxFlux columns from target model file
-            minFlux = data.MinFlux;
-            maxFlux = data.MaxFlux;
+            minFluxTarget = data.MinFlux;
+            maxFluxTarget = data.MaxFlux;
 
             % Create and fill cell arrays (min, max) with rounded flux
             % values from target model
-            minFlux_r = cell(height(minFlux),1);
-            maxFlux_r = cell(height(maxFlux),1);
+            minFluxTarget_r = cell(height(minFluxTarget),1);
+            maxFluxTarget_r = cell(height(maxFluxTarget),1);
 
             for k=1:1:height(data)
-                minFlux_r{k} = round(str2double(minFlux{k}),6);
-                maxFlux_r{k} = round(str2double(maxFlux{k}),6);
+                minFluxTarget_r{k} = round(str2double(minFluxTarget{k}),6);
+                maxFluxTarget_r{k} = round(str2double(maxFluxTarget{k}),6);
             end
 
             % Create and fill flux ratio cell arrays (target/source)
-            MinFluxRatio = cell(height(minFlux),1);
-            MaxFluxRatio = cell(height(minFlux),1);
+            MinFluxRatio = cell(height(minFluxTarget),1);
+            MaxFluxRatio = cell(height(minFluxTarget),1);
 
             % Calculate ratios for all reaction fluxes
-            for n=1:1:length(minFlux)
-                minRatio = minFlux_r{n}/minFluxSource_r{n};
-                maxRatio = maxFlux_r{n}/maxFluxSource_r{n};
+            for n=1:1:length(minFluxTarget)
+                minRatio = minFluxTarget_r{n}/minFluxSource_r{n};
+                maxRatio = maxFluxTarget_r{n}/maxFluxSource_r{n};
                 minRatio = round(minRatio,6);
                 maxRatio = round(maxRatio,6);
 
                 % Source flux is 0
-                if minFluxSource_r{n} == 0 && minFlux_r{n} ~= 0
+                if minFluxSource_r{n} == 0 && minFluxTarget_r{n} ~= 0
                     minRatio = 'up';
                 end
-                if maxFluxSource_r{n} == 0 && maxFlux_r{n} ~= 0
+                if maxFluxSource_r{n} == 0 && maxFluxTarget_r{n} ~= 0
                     maxRatio = 'up';
                 end
 
                 % Target flux is 0
-                if minFlux_r{n} == 0 &&  minFluxSource_r{n} ~= 0
+                if minFluxTarget_r{n} == 0 &&  minFluxSource_r{n} ~= 0
                     minRatio = 'down';
                 end
-                if maxFlux_r{n} == 0 && maxFluxSource_r{n} ~= 0
+                if maxFluxTarget_r{n} == 0 && maxFluxSource_r{n} ~= 0
                     maxRatio = 'down';
                 end
 
                 % Both are 0
-                if minFlux_r{n} == 0 &&  minFluxSource_r{n} == 0
+                if minFluxTarget_r{n} == 0 &&  minFluxSource_r{n} == 0
                     minRatio = 1;
                 end
-                if maxFlux_r{n} == 0 && maxFluxSource_r{n} == 0
+                if maxFluxTarget_r{n} == 0 && maxFluxSource_r{n} == 0
                     maxRatio = 1;
                 end
 
@@ -108,7 +109,7 @@ function calculateFluxShifts(source, target)
             end
 
             % Save results
-            folderName = 'resultsPostOptimization\fluxShifts\';
+            folderName = 'resultsPostOptimization/fluxShifts/';
             if ~exist(folderName, 'dir')
                mkdir(folderName)
             end
@@ -117,7 +118,7 @@ function calculateFluxShifts(source, target)
             ratioData.Properties.VariableNames = {strcat('MinFlux_',source{1}), strcat('MaxFlux_',source{1}), 'MinFluxRatio', 'MaxFluxRatio'};
             newData = [data ratioData];
 
-            temp = split(T(i).name,".");
+            temp = split(xlsFilesInfo(i).name,".");
             resultFileName = strcat(temp(1),'_flux_shifts.xls');
             fullResultFilePath = strcat(folderName, resultFileName);
             writetable(newData,fullResultFilePath{1},'AutoFitWidth',false);

src/analysis/exploration/biomassPrecursorCheck.m

@@ -32,6 +32,7 @@
 %                    (returned only if checkConservedQuantities = true)
 %
 % .. Authors: - Pep Charusanti & Richard Que (July 2010)
+%             - Hadjar Rahou (Nov 2023) (fix a bug)
 
 if ~exist('checkCoupling','var') || isempty(checkCoupling)
     checkCoupling = 0;
@@ -118,12 +119,19 @@
         EMV = findElementaryMoietyVectors(model);
     else
         % atom transition network is supplied.
-        EMV = identifyConservedMoieties(model, ATN);
-        types = classifyMoieties(EMV, model.S);
-        EMV = EMV(:, strcmp(types, 'Internal'));
+        %EMV = identifyConservedMoieties(model, ATN);
+        %types = classifyMoieties(EMV, model.S);
+        %EMV = EMV(:, strcmp(types, 'Internal'));
+        arm = identifyConservedMoieties(model, ATN); %Hadjar
+        L=arm.L; 
+        types = classifyMoieties(L, model.S);
+        L = L(strcmp(types, 'Internal'),:);
+        EMV=L';
     end
     % biomass metabolites that contain conserved moieties
-    mCofactor = any(model.S(:, colS_biomass) ~= 0, 2) & any(EMV, 2);
+    %mCofactor = any(model.S(:, colS_biomass) ~= 0, 2) & any(EMV, 2);
+    mCofactor = (any(model.S(:, colS_biomass) < 0, 2) & any(EMV, 2)); %Hadjar
+    %pause(2);
     % elementary moieties involved in biomass production
     cofactorPairMatrix = EMV(:, any(EMV(mCofactor, :), 1));
     % each cell of cofactorPair stores the set of cofactor metabolites to be produced.

src/analysis/exploration/printConstraints.m

@@ -18,7 +18,7 @@ function printConstraints(model, minInf, maxInf, rxnSelection, modelAfter, print
 %
 
 % .. Authors:
-%       - Ines Thiele 02/09, Ronan Fleming 2020
+%       - Ines Thiele 02/09, Ronan Fleming 2020, Yanjun Liu 2023
 
 if ~exist('minInf','var')|| isempty(minInf)
     minInf=-Inf;
@@ -67,7 +67,7 @@ function printConstraints(model, minInf, maxInf, rxnSelection, modelAfter, print
 
 
 if ~any(closedRxnBool | reversibleRxnBool | fwdRxnBool0b | revRxnBool0b | fwdRxnBoolNon0b | revRxnBoolNon0b)
-    boolRemainder = rxnSelection & ~(closedRxnBool | reversibleRxnBool | fwdRxnBool | revRxnBool);
+    boolRemainder = rxnSelection & ~(closedRxnBool | reversibleRxnBool | fwdRxnBool0b | revRxnBool0b | fwdRxnBoolNon0b | revRxnBoolNon0b);
     warning('no subset with bounds between [minInf maxInf]')
 else
     boolRemainder=0;

src/analysis/exploration/printFluxVector.m

@@ -20,7 +20,8 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
 %
 
 % .. Authors:
-%       - Markus Herrgard, Ronan Fleming
+%    - Markus Herrgard, Ronan Fleming
+%    - Farid Zare    21/11/2023    Corrected outputing multiple rxn names
 
 if isempty(fluxData)
     return
@@ -70,7 +71,6 @@ function printFluxVector(model, fluxData, nonZeroFlag, excFlag, sortCol, fileNam
 % Add reaction formulas
 if formulaFlag
     if nonZeroFlag
-        labels = [model.rxns,model.rxns];
         %only generate the formulas for the nonzero entries
         labels(bool,end+1) = printRxnFormula(model, labels(bool), false, false);
     else

src/analysis/filterNonFluxReactions/filterNonFluxReactions.m

@@ -16,8 +16,9 @@ function filterNonFluxReactions(phenotype)
 % .. Authors:
 %       - Kristina Grausa 05/16/2022
 %       - Kristina Grausa 08/23/2022 - standard header and formatting
+%       - Farid Zare      20/11/2023 - Repository addresses are correcteds
 
-    dest = string(strcat('resultsPostOptimization\contextSpecificModels\*', phenotype, '*.*'));
+    dest = string(strcat('resultsPostOptimization/contextSpecificModels/*', phenotype, '*.*'));
     S = dir(dest);
     count = length(S);
 
@@ -27,7 +28,7 @@ function filterNonFluxReactions(phenotype)
 
     for i=1:1:count
         sheetName = 'Reaction List';
-        filename = strcat(S(i).folder, '\', S(i).name);
+        filename = strcat(S(i).folder, '/', S(i).name);
         data=readtable(filename,'Sheet',sheetName);
 
         % Copy headers
@@ -75,7 +76,7 @@ function filterNonFluxReactions(phenotype)
     end
 
     % Save results
-    folderName = 'resultsPostOptimization\nonFluxReactions\';
+    folderName = 'resultsPostOptimization/nonFluxReactions/';
     if ~exist(folderName, 'dir')
        mkdir(folderName)
     end

src/analysis/filterRateLimittingReactions/filterRateLimittingReactions.m

@@ -7,17 +7,19 @@ function filterRateLimittingReactions(phenotype)
 %
 % USAGE:
 %
-%   filterRateLimittingReactions(phenotype)
+%    filterRateLimittingReactions(phenotype)
 %
 % INPUTS:
-%   phenotype:              char representing the phenotype name provided
-%                           in each model name of the same phanotype
+%    phenotype:              char representing the phenotype name provided
+%                            in each model name of the same phanotype
 %
 % .. Authors:
 %       - Kristina Grausa 05/16/2022
 %       - Kristina Grausa 08/23/2022 - standard header and formatting
+%       - Farid Zare      11/21/2023 - Correction of repository address format 
+%
 
-    dest = string(strcat('resultsPostOptimization\contextSpecificModels\*', phenotype, '*.xls'));
+    dest = string(strcat('resultsPostOptimization/contextSpecificModels/*', phenotype, '*.xls'));
     S = dir(dest);
     count = length(S);
 
@@ -27,7 +29,7 @@ function filterRateLimittingReactions(phenotype)
 
     for i=1:1:count
         sheetName = 'Reaction List';
-        filename = strcat(S(i).folder, '\', S(i).name);
+        filename = strcat(S(i).folder, '/', S(i).name);
         data=readtable(filename,'Sheet',sheetName);
 
         % Copy headers
@@ -79,7 +81,7 @@ function filterRateLimittingReactions(phenotype)
     end
 
     % Save results
-    folderName = 'resultsPostOptimization\rateLimittingReactions\';
+    folderName = 'resultsPostOptimization/rateLimittingReactions/';
     if ~exist(folderName, 'dir')
        mkdir(folderName)
     end

src/analysis/multiSpecies/SteadyCom/subroutines/SteadyComSubroutines.m

@@ -234,7 +234,7 @@
         spName = spAbbr;
     end
     if isfield(infoCom, 'spBm')
-        indCom.spBm = modelCom.rxns(infoCom.spBm);
+        indCom.spBm = modelCom.rxns(infoCom.spBm); %modelCom.infoCom.spBm
     end
     if isfield(infoCom,'spATPM')
         indCom.spATPM = modelCom.rxns(infoCom.spATPM);

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

@@ -22,10 +22,13 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
 %                       statistical analysis are saved
 % violinPath:           char with path of directory where violin plots are
 %                       saved
+% plotType              Type of plot to be created: violin plot (default)
+%                       or boxplot. Allowed entries: "ViolinPlot","Boxplot"
 % sampleGroupHeaders    list of one or more column headers in file with the
 %                       sample information that should be analyzed 
 %                       (e.g., disease status, age). If not provided, the
 %                       second column will be used.
+% 
 %
 % .. Author: Almut Heinken, 12/2020
 
@@ -34,6 +37,7 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
 parser.addRequired('resPath', @ischar);
 parser.addParameter('statPath', [pwd filesep 'Statistics'], @ischar);
 parser.addParameter('violinPath', [pwd filesep 'ViolinPlots'], @ischar);
+parser.addParameter('plotType', 'ViolinPlot', @ischar);
 parser.addParameter('sampleGroupHeaders', '', @iscellstr);
 
 parser.parse(infoFilePath, resPath, varargin{:});
@@ -42,8 +46,13 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
 resPath = parser.Results.resPath;
 statPath = parser.Results.statPath;
 violinPath = parser.Results.violinPath;
+plotType = parser.Results.plotType;
 sampleGroupHeaders = parser.Results.sampleGroupHeaders;
 
+if ~any(strcmp(plotType,{'ViolinPlot','Boxplot'}))
+    error('Invalid entry for plot type!')
+end
+
 % create the folders
 mkdir(statPath)
 mkdir(violinPath)
@@ -108,13 +117,13 @@ function analyzeMgPipeResults(infoFilePath,resPath,varargin)
                 writetable(cell2table(significantFeatures),[statPath filesep filename '_' sampleGroupHeaders{j} '_SignificantFeatures'],'FileType','text','WriteVariableNames',false,'Delimiter','tab');
             end
 
-            % create violin plots
+            % create plots
             currentDir=pwd;
             cd(violinPath)
 
             % create violin plots for net uptake and secretion files
             if any(contains(fileList{i,1},{'net_uptake_fluxes.csv','net_secretion_fluxes.csv'}))
-                makeViolinPlots(sampleData, infoFile, 'stratification',sampleGroupHeaders{j}, 'plottedFeature', filename, 'unit', 'mmol/person/day')
+                makeViolinPlots(sampleData, infoFile, 'stratification',sampleGroupHeaders{j}, 'plotType' , plotType, 'plottedFeature', filename, 'unit', 'mmol/person/day')
             end
             cd(currentDir)
         end