changeset 295:626b6d1de075 draft

Uploaded

--- a/COBRAxy/flux_to_map.py	Thu May 15 18:30:02 2025 +0000
+++ b/COBRAxy/flux_to_map.py	Fri May 16 10:56:01 2025 +0000
@@ -65,6 +65,11 @@
         type = float, 
         default = 0.1, 
         help = 'P-Value threshold (default: %(default)s)')
+
+    parser.add_argument(
+        '-adj' ,'--adjusted',
+        type = utils.Bool("adjusted"), default = False, 
+        help = 'Apply the FDR (Benjamini-Hochberg) correction (default: %(default)s)')
     
     parser.add_argument(
         '-fc', '--fChange',
@@ -72,7 +77,6 @@
         default = 1.5, 
         help = 'Fold-Change threshold (default: %(default)s)')
     
-
     parser.add_argument(
         '-op', '--option',
         type = str, 
@@ -649,7 +653,7 @@
 
 def compareDatasetPair(dataset1Data :List[List[float]], dataset2Data :List[List[float]], ids :List[str]) -> Tuple[Dict[str, List[Union[float, FoldChange]]], float]:
     #TODO: the following code still suffers from "dumbvarnames-osis"
-    tmp :Dict[str, List[Union[float, FoldChange]]] = {}
+    comparisonResult :Dict[str, List[Union[float, FoldChange]]] = {}
     count   = 0
     max_z_score = 0
     for l1, l2 in zip(dataset1Data, dataset2Data):
@@ -664,10 +668,22 @@
             f_c = fold_change(avg1, avg2)
             if np.isfinite(z_score) and max_z_score < abs(z_score): max_z_score = abs(z_score)
             
-            tmp[reactId] = [float(p_value), f_c, z_score, avg1, avg2]
+            comparisonResult[reactId] = [float(p_value), f_c, z_score, avg1, avg2]
         except (TypeError, ZeroDivisionError): continue
+
+    # Apply multiple testing correction if set by the user
+    if ARGS.adjusted:
+        
+        # Retrive the p-values from the comparisonResult dictionary
+        reactIds = list(comparisonResult.keys())
+        pValues = [comparisonResult[reactId][0] for reactId in reactIds]
+        
+        # Apply the Benjamini-Hochberg correction and update
+        adjustedPValues = st.multipletests(pValues, method="fdr_bh")[1]
+        for i, reactId in enumerate(reactIds):
+            comparisonResult[reactId][0] = adjustedPValues[i]
     
-    return tmp, max_z_score
+    return comparisonResult, max_z_score
 
 def computeEnrichment(class_pat :Dict[str, List[List[float]]], ids :List[str]) -> List[Tuple[str, str, dict, float]]:
     """