Changed how flux reduction is calculated

scripts/knockout_fva.py

@@ -38,7 +38,7 @@ def remove_blocked_reactions(model, proc, progress):
         print("No blocked reactions could be found.")
 
 
-def run_fva(model, case, rxns, eps, rxnsOfInterest, proc, fraction_opt, reduce):
+def run_fva(model, case, rxns, eps, rxnsOfInterest, proc, fraction_opt, reduce, WT):
     if case == 'forward':
         with model:
             for r in rxns:  # Set all the reaction bounds for the forward FVA
@@ -93,8 +93,28 @@ def run_fva(model, case, rxns, eps, rxnsOfInterest, proc, fraction_opt, reduce):
                                                    processes=proc)
             else:
                 for r in rxns:  # Reduce all the reaction bounds for the mutant forward FVA
+                    # calculate new bounds
+                    # Check sign of WT bounds:
                     rxn = model.reactions.get_by_id(r)
-                    model.reactions.get_by_id(rxn.id).upper_bound = model.reactions.get_by_id(rxn.id).lower_bound = reduce*eps
+                    pd.set_option('display.max_rows', 1000)
+                    # print(WT)
+                    old_ub = WT.loc[rxn.id][1]
+                    old_lb = WT.loc[rxn.id][0]
+                    if old_ub < 0:  # Ex: [-20, -5]
+                        new_lb = old_lb - ((old_ub - old_lb)*reduce)  # -20 - ((-20+5)*0.5) = - 12.5
+                        new_ub = old_ub  # No difference
+                    elif old_lb > 0:  # Ex: [100, 500]
+                        new_lb = old_lb  # No difference
+                        new_ub = old_ub - ((old_ub - old_lb) * reduce)  # 500 - ((500-100)*0.5) = 300
+
+                    elif old_lb == 0 and old_ub == 0:
+                        new_ub = 0
+                        new_lb = 0
+                    else:  # Ex: [-50, 10]
+                        new_lb = old_lb - ((old_ub - abs(old_lb)) * reduce)/2  # -50 - ((10-50)*0.5)/2 = -40
+                        new_ub = old_ub - ((old_ub - abs(old_lb)) * reduce)/2  # 10 - ((10-50)*0.5)/2 = 0
+                    model.reactions.get_by_id(rxn.id).upper_bound = new_ub
+                    model.reactions.get_by_id(rxn.id).lower_bound = new_lb
                 try:
                     mutantf = flux_variability_analysis(model, reaction_list=rxnsOfInterest,
                                                         fraction_of_optimum=fraction_opt, processes=proc)