Add files via upload

BigDataInSilicoBiologyGroup · Apr 11, 2022 · b4a2161 · b4a2161
1 parent df6ccf0
commit b4a2161
Show file tree

Hide file tree

Showing 17 changed files with 170 additions and 95 deletions.
diff --git a/Scripts/TestCase_calculateMinimumRequirements.m b/Scripts/TestCase_calculateMinimumRequirements.m
@@ -1,13 +1,14 @@
 modelPath = 'Data/Yeast_8_4_0.xls';
 trDataPath = 'Data/TranscriptomicsData.xlsx';
 mediumDataPath = 'Data/MediumData.xlsx';
-growthNotAffectingGeneDel = 1;
+growthNotAffectingGeneDel = 0;
 thApproach = 3;
 lowerTh = 5;
 upperTh = 50;
 objective = 'r_2111';
+percentile = 0;
 try
-    calculateMinimumRequirements(modelPath, trDataPath, mediumDataPath, growthNotAffectingGeneDel, thApproach, lowerTh, upperTh, objective);
+    calculateMinimumRequirements(modelPath, trDataPath, mediumDataPath, growthNotAffectingGeneDel, thApproach, lowerTh, upperTh, objective, percentile);
     disp('calculateMinimumRequirements.m... OK');
 catch e
     disp('calculateMinimumRequirements.m... NOK');

diff --git a/Scripts/TestCase_createContextSpecificModel.m b/Scripts/TestCase_createContextSpecificModel.m
@@ -1,17 +1,19 @@
 modelPath = 'Data/Yeast_8_4_0.xls';
-trDataPath = 'Data/TranscriptomicData.xlsx';
+trDataPath = 'Data/TranscriptomicsData.xlsx';
 mediumDataPath = 'Data/MediumData.xlsx';
 growthNotAffectingGeneDel = 1;
 meetMinimumReq = 0;
 thApproach = 3;
 lowerTh = 2;
 upperTh = 25;
 objective = 'r_2111';
-gmMAX = 1;
+gmAndOperation = 'MIN';
+gmOrOperation = 'MAX';
 constrAll = 0;
 excludeBiomassEq = 0;
+percentile=0;
 try
-    model = createContextSpecificModel(modelPath, trDataPath, mediumDataPath, growthNotAffectingGeneDel, meetMinimumReq, thApproach, lowerTh, upperTh, objective, gmMAX, constrAll, excludeBiomassEq);
+    model = createContextSpecificModel(modelPath, trDataPath, mediumDataPath, growthNotAffectingGeneDel, meetMinimumReq, thApproach, lowerTh, upperTh, objective, gmAndOperation, gmOrOperation, constrAll, excludeBiomassEq,percentile);
     disp('createContextSpecificModel.m... OK');
 catch e
     disp('createContextSpecificModel.m... NOK');

diff --git a/Scripts/TestCase_deleteInactiveGenes.m b/Scripts/TestCase_deleteInactiveGenes.m
@@ -5,9 +5,10 @@
 lowerTh=0;
 upperTh=30;
 sheetIndex=1;
-growthNotAffectingGeneDel=1;
+growthNotAffectingGeneDel=0;
+percentile=0;
 try
-    model = deleteInactiveGenes(model, trData, trDataPath, thApproach, lowerTh, upperTh, sheetIndex, growthNotAffectingGeneDel);
+    model = deleteInactiveGenes(model, trData, trDataPath, thApproach, lowerTh, upperTh, sheetIndex, growthNotAffectingGeneDel,percentile);
     disp('deleteInactiveGenes.m... OK');
 catch e 
     disp('deleteInactiveGenes.m... NOK');

diff --git a/Scripts/TestCase_fluxShifts.m b/Scripts/TestCase_fluxShifts.m
@@ -8,12 +8,12 @@
 threshold = 0; % (possible: 0-5)
 
 % Find the files described in the comments above
-dest = strcat('..\Results post-optimization\*S47D*.xls');
+dest = strcat('..\Results post-optimization\Context-specific models\*S47D*.xls');
 S = dir(dest);
 count = length(S);
 
 % Find WT data file (is used to compare with mutants)
-destWT = strcat('..\Results post-optimization\SRR8994357_WT_GT1_0.xls');
+destWT = strcat('..\Results post-optimization\Context-specific models\SRR8994357_WT_GT1_0.xls');
 WT = dir(destWT);
 sheetName = 'Reaction List';
 filename = strcat(WT.folder, '\', WT.name);

diff --git a/Scripts/calculateFluxShifts.m b/Scripts/calculateFluxShifts.m
@@ -36,6 +36,7 @@ function calculateFluxShifts(source, target)
             filename = strcat(T(i).folder, '\', T(i).name); % full path
             data=readtable(filename,'Sheet',sheetName);
 
+            try
             % Read target fluxes
             minFlux = data.MinFlux;
             maxFlux = data.MaxFlux;
@@ -102,6 +103,9 @@ function calculateFluxShifts(source, target)
             resultFileName = strcat(temp(1),'_flux_shifts.xls');
             fullResultFilePath = strcat(folderName, resultFileName);
             writetable(newData,fullResultFilePath{1},'AutoFitWidth',false);
+            catch e
+                disp(e);
+            end
         end
     end
 end
diff --git a/Scripts/calculateMinimumRequirements.m b/Scripts/calculateMinimumRequirements.m
@@ -6,17 +6,20 @@ function calculateMinimumRequirements(modelPath, trDataPath, mediumDataPath, gro
         model = modelOriginal;
         trData=readtable(trDataPath,'Sheet',trSheets{i}); 
         % Calculate percentile
-        if percentile == 1        
-            lowerTh = calculatePercentile(trData.Data, lowerTh); % TEST
-            upperTh = calculatePercentile(trData.Data, upperTh);
+        if percentile == 1      
+            lowerThreshold = calculatePercentile(trData.Data, lowerTh);
+            upperThreshold = calculatePercentile(trData.Data, upperTh);
+        else
+            lowerThreshold = lowerTh;
+            upperThreshold = upperTh;
         end
         inactiveGenes = {};
         if thApproach == 1
-            inactiveGenes = findGenesBelowThresholdGT1(lowerTh, trData.Geneid, trData.Data);
+            inactiveGenes = findGenesBelowThresholdGT1(lowerThreshold, trData.Geneid, trData.Data);
         elseif thApproach == 2
-            inactiveGenes = findGenesBelowThresholdLocal1(lowerTh,trDataPath,i);
+            inactiveGenes = findGenesBelowThresholdLocal1(lowerThreshold,trDataPath,i);
         else
-            inactiveGenes = findGenesBelowThresholdLocal2(lowerTh, upperTh, trDataPath,i);
+            inactiveGenes = findGenesBelowThresholdLocal2(lowerThreshold, upperThreshold, trDataPath,i);
         end
 
         % Non expressed genes
@@ -29,6 +32,7 @@ function calculateMinimumRequirements(modelPath, trDataPath, mediumDataPath, gro
                         try
                             % Test is gene deletions affects growth
                             [grRatio, grRateKO, grRateWT, hasEffect, delRxns] = singleGeneDeletion(model, 'FBA', inactiveGenes(j));
+                            disp(inactiveGenes(j));
                             if grRatio == 1
                                 genes_to_delete{counter} = inactiveGenes{j};
                                 counter = counter + 1;

diff --git a/Scripts/calculatePercentile.m b/Scripts/calculatePercentile.m
@@ -1,14 +1,14 @@
 function value = calculatePercentile(expressionValues, k)
-    if ~isempty(k)
+    if ~isnan(k)
         C = unique(expressionValues,'sorted');
-        n = length(C); % try with original length as well
-        Lp =  (n+1) * (str2double(k)/100); % Position Locator
+        n = length(C);
+        Lp =  (n+1) * (k/100); % Position Locator
         if round(Lp) == Lp
             value = C(Lp);
         else
             value = round(C(floor(Lp)) + (Lp-floor(Lp)) * (C(ceil(Lp))-C(floor(Lp))));
         end
     else 
-        value = k; % empty
+        value = k;
     end
 end
diff --git a/Scripts/createContextSpecificModel.m b/Scripts/createContextSpecificModel.m
@@ -1,11 +1,10 @@
-function model = createContextSpecificModel(modelPath, trDataPath, mediumDataPath, growthNotAffectingGeneDel, meetMinimumReq, thApproach, lowerTh, upperTh, objective, gmMAX, constrAll, excludeBiomassEq, biomassId)
+function model = createContextSpecificModel(modelPath, trDataPath, mediumDataPath, growthNotAffectingGeneDel, meetMinimumReq, thApproach, lowerTh, upperTh, objective, gmAndOperation, gmOrOperation, constrAll, excludeBiomassEq, biomassId, percentile)
     % model - ...
     % geneNames - ...
     % ExpressionValues - ...
     % geneMapping - 1 for AND/OR=MIN/MAX; 2 for AND/OR=MIN/SUM
     % constrAll - true for constrain all reactions; 2 for constrain only
-    % irreversible reactions
-
+    % irreversible reactions   
 
     trSheets = sheetnames(trDataPath);
     modelOriginal=readCbModel(modelPath);
@@ -25,7 +24,7 @@
             try
                 minRequirements = readtable(strcat('Results post-optimization/Minimum requirements/',trSheets{s},'.xls')).FBAMin; 
             catch e
-                minRequirements = array2table(zeros(length(model.rxns),1));
+                minRequirements = {};
             end
         end
 
@@ -68,7 +67,11 @@
                             andsToCompare = [andsToCompare tr_Val];
                         end
                     end
-                    orsToCompare = [orsToCompare min(andsToCompare)];
+                    if strcmp(gmAndOperation,'MIN')
+                        orsToCompare = [orsToCompare min(andsToCompare)];
+                    else 
+                        orsToCompare = [orsToCompare round(geomean(andsToCompare))];
+                    end
                 end
 
                 operands = split(model.grRules{i});
@@ -79,8 +82,7 @@
                     end
                 end
                 log{i,4} = value;
-
-                if gmMAX == 1
+                if strcmp(gmOrOperation,'MAX')
                     log{i,5} = max(orsToCompare);
                 else 
                     log{i,5} = sum(orsToCompare);
@@ -101,7 +103,11 @@
                     end
                 end
                 log{i,4} = value;
-                log{i,5} = min(numbersToCompare);
+                if strcmp(gmAndOperation,'MIN')
+                    log{i,5} = min(numbersToCompare);
+                else 
+                    log{i,5} = round(geomean(numbersToCompare));
+                end
 
             % Gene combination OR
             elseif contains(model.grRules{i},"or")
@@ -118,7 +124,7 @@
                     end
                 end
                 log{i,4} = value;
-                if gmMAX == 1
+                if strcmp(gmOrOperation,'MAX')
                     log{i,5} = max(orsToCompare);
                 else 
                     log{i,5} = sum(orsToCompare);
@@ -127,21 +133,32 @@
 
             %Set reaction bounds
             if ~isempty(log{i,5})
-                if (meetMinimumReq == 1 && log{i,5} >= str2double(minRequirements{i})) || meetMinimumReq == 0
-                    if constrAll == 1 && model.lb(i) < 0
-                        value = log{i,5};
-                        model = changeRxnBounds(model,model.rxns{i},-value,'l');
-                        model = changeRxnBounds(model,model.rxns{i},value,'u');
-                    elseif model.lb(i) == 0
-                        model = changeRxnBounds(model,model.rxns{i},log{i,5},'u');
+                try
+                    constrainRxn = true;
+                    if ~isempty(minRequirements)
+                        if log{i,5} < str2double(minRequirements{i})
+                            constrainRxn = false;
+                        end
+                    end
+
+                    if constrainRxn
+                        if constrAll == 1 && model.lb(i) < 0
+                            value = log{i,5};
+                            model = changeRxnBounds(model,model.rxns{i},-value,'l');
+                            model = changeRxnBounds(model,model.rxns{i},value,'u');
+                        elseif model.lb(i) == 0
+                            model = changeRxnBounds(model,model.rxns{i},log{i,5},'u');
+                        end
                     end
+                catch e
+                  disp(e);
                 end
             end
         end
 
 
         % Delete Inactive genes
-        model = deleteInactiveGenes(model, trData, trDataPath, thApproach, lowerTh, upperTh, s, growthNotAffectingGeneDel);
+        model = deleteInactiveGenes(model, trData, trDataPath, thApproach, lowerTh, upperTh, s, growthNotAffectingGeneDel, percentile);
 
 
         % Apply medium data

diff --git a/Scripts/deleteInactiveGenes.m b/Scripts/deleteInactiveGenes.m
@@ -1,5 +1,9 @@
-function model = deleteInactiveGenes(model, trData, trDataPath, thApproach, lowerTh, upperTh, sheetIndex, growthNotAffectingGeneDel)      
-
+function model = deleteInactiveGenes(model, trData, trDataPath, thApproach, lowerTh, upperTh, sheetIndex, growthNotAffectingGeneDel, percentile)      
+    % Calculate percentile
+    if percentile == 1        
+        lowerTh = calculatePercentile(trData.Data, lowerTh);
+        upperTh = calculatePercentile(trData.Data, upperTh);
+    end
     inactiveGenes = {};
     if thApproach == 1
         inactiveGenes = findGenesBelowThresholdGT1(lowerTh, trData.Geneid, trData.Data);
@@ -18,6 +22,7 @@
                 if growthNotAffectingGeneDel == 1
                     % Test is gene deletions affects growth
                     [grRatio, grRateKO, grRateWT, hasEffect, delRxns] = singleGeneDeletion(model, 'FBA', inactiveGenes(j));
+                    disp(inactiveGenes(j));
                     if grRatio == 1
                         genes_to_delete{counter} = inactiveGenes{j};
                         counter = counter + 1;

diff --git a/Scripts/filterRateLimittingReactions.m b/Scripts/filterRateLimittingReactions.m
@@ -17,12 +17,16 @@ function filterRateLimittingReactions(phenotype)
         counter = 2;
         for n=1:1:height(data)
             ub = str2double(data.UpperBound{n});
+            try
             maxFlux = str2double(data.MaxFlux{n});
             gpr = data.GPR{n};
             if ub == maxFlux && ub ~= 0 && maxFlux ~= 1000 && ~isempty(gpr)
                 result(counter,:) = table2cell(data(n,:));
                 counter = counter + 1;
             end
+            catch e
+             disp(e);
+            end
         end
         resultAll{i} = result;
         result = {};

diff --git a/Scripts/findGenesAboveThresholdLocal1.m b/Scripts/findGenesAboveThresholdLocal1.m
@@ -11,8 +11,8 @@
         % Get target trancriptomics sheet
         trTarget=readtable(trDataPath,'Sheet',trSheets{sheetIndex}); 
 
-        % Get all data sets of the same phemotype    
-        for i=1:1:height(trSheets) % GENES MUST BE IN THE SAME ORDER 
+        % Get all data sets of the same phemotype (GENES MUST BE IN THE SAME ORDER)  
+        for i=1:1:height(trSheets) 
             if contains(trSheets{i}, phenotype)
                 data=readtable(trDataPath,'Sheet',trSheets{i}); 
                 if length(trDataAll) == 0

diff --git a/Scripts/findGenesAboveThresholdLocal2.m b/Scripts/findGenesAboveThresholdLocal2.m
@@ -11,8 +11,9 @@
         % Get target trancriptomics sheet
         trTarget=readtable(trDataPath,'Sheet',trSheets{sheetIndex}); 
 
-        % Get all data sets of the same phemotype    
-        for i=1:1:height(trSheets) % GENES MUST BE IN THE SAME ORDER 
+        % Get all data sets of the same phemotype (GENES MUST BE IN THE
+        % SAME ORDER)
+        for i=1:1:height(trSheets)
             if contains(trSheets{i}, phenotype)
                 data=readtable(trDataPath,'Sheet',trSheets{i}); 
                 if length(trDataAll) == 0

diff --git a/Scripts/findGenesBelowThresholdLocal1.m b/Scripts/findGenesBelowThresholdLocal1.m
@@ -3,16 +3,16 @@
         trDataAll = {};
         cnt = 1;
 
-        % Get transcriptomics data sheet names
+        % Get transcriptomics data sheet (phenotype) name
         trSheets = sheetnames(trDataPath);
         sheetNameParts = split(trSheets{sheetIndex},'_');
         phenotype = sheetNameParts{2};
 
-        % Get target trancriptomics sheet
+        % Get target trancriptomics sheet data
         trTarget=readtable(trDataPath,'Sheet',trSheets{sheetIndex}); 
 
-        % Get all data sets of the same phemotype    
-        for i=1:1:height(trSheets) % GENES MUST BE IN THE SAME ORDER 
+        % Get all data sets of the same phemotype (GENES MUST BE IN THE SAME ORDER)  
+        for i=1:1:height(trSheets) 
             if contains(trSheets{i}, phenotype)
                 data=readtable(trDataPath,'Sheet',trSheets{i}); 
                 if length(trDataAll) == 0

diff --git a/Scripts/findGenesBelowThresholdLocal2.m b/Scripts/findGenesBelowThresholdLocal2.m
@@ -11,8 +11,8 @@
         % Get target trancriptomics sheet
         trTarget=readtable(trDataPath,'Sheet',trSheets{sheetIndex}); 
 
-        % Get all data sets of the same phemotype    
-        for i=1:1:height(trSheets) % GENES MUST BE IN THE SAME ORDER 
+        % Get all data sets of the same phemotype (GENES MUST BE IN THE SAME ORDER)    
+        for i=1:1:height(trSheets)
             if contains(trSheets{i}, phenotype)
                 data=readtable(trDataPath,'Sheet',trSheets{i}); 
                 if length(trDataAll) == 0

diff --git a/Scripts/findTranscriptionValue.m b/Scripts/findTranscriptionValue.m
@@ -1,9 +1,6 @@
 function value = findTranscriptionValue(geneId, trData)
     geneIds = trData.Geneid;
     expressionValues = trData.Data;
-    for k=1:1:length(geneIds)
-       if strcmp(geneId, geneIds(k))
-           value = expressionValues(k);
-       end
-   end
+    index = find(strcmp(trData.Geneid, geneId));
+    value = expressionValues(index);
 end