diff COBRAxy/flux_to_map.py @ 302:753347af8bc9 draft

Uploaded

--- a/COBRAxy/flux_to_map.py	Tue May 20 14:47:56 2025 +0000
+++ b/COBRAxy/flux_to_map.py	Tue May 20 14:49:15 2025 +0000
@@ -271,8 +271,8 @@
     """
     Encodes possible arrow colors based on their meaning in the enrichment process.
     """
-    Invalid       = "#BEBEBE" # gray, fold-change under treshold
-    Transparent   = "#ffffff00" # white, not significant p-value
+    Invalid       = "#BEBEBE" # gray, fold-change under treshold or not significant p-value
+    Transparent   = "#ffffff00" # transparent, to make some arrow segments disappear
     UpRegulated   = "#ecac68" # red, up-regulated reaction
     DownRegulated = "#6495ed" # blue, down-regulated reaction
 
@@ -625,12 +625,18 @@
             # Perform Kolmogorov-Smirnov test
             _, p_value = st.ks_2samp(dataset1Data, dataset2Data)
         case "ttest_p":
+            # Datasets should have same size
+            if len(dataset1Data) != len(dataset2Data):
+                raise ValueError("Datasets must have the same size for paired t-test.")
             # Perform t-test for paired samples
             _, p_value = st.ttest_rel(dataset1Data, dataset2Data)
         case "ttest_ind":
             # Perform t-test for independent samples
             _, p_value = st.ttest_ind(dataset1Data, dataset2Data)
         case "wilcoxon":
+            # Datasets should have same size
+            if len(dataset1Data) != len(dataset2Data):
+                raise ValueError("Datasets must have the same size for Wilcoxon signed-rank test.")
             # Perform Wilcoxon signed-rank test
             _, p_value = st.wilcoxon(dataset1Data, dataset2Data)
         case "mw":
@@ -674,14 +680,19 @@
     # 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]
+        # Retrieve the p-values from the comparisonResult dictionary, they have to be different from NaN
+        validPValues = [(reactId, result[0]) for reactId, result in comparisonResult.items() if not np.isnan(result[0])]
+
+        if not validPValues:
+            return comparisonResult, max_z_score
         
-        # 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]
+        # Unpack the valid p-values
+        reactIds, pValues = zip(*validPValues)
+        # Adjust the p-values using the Benjamini-Hochberg method
+        adjustedPValues = st.false_discovery_control(pValues)
+        # Update the comparisonResult dictionary with the adjusted p-values
+        for reactId , adjustedPValue in zip(reactIds, adjustedPValues):
+            comparisonResult[reactId][0] = adjustedPValue
     
     return comparisonResult, max_z_score