From 8f80e02f389e9211ea5f7dfdc00978c997b650ba Mon Sep 17 00:00:00 2001
From: haydensun <hayden_sun@hotmail.com>
Date: Fri, 15 May 2020 00:08:53 +0800
Subject: [PATCH] Add stats module

---
 mamotif/stats.py | 57 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 57 insertions(+)
 create mode 100644 mamotif/stats.py

diff --git a/mamotif/stats.py b/mamotif/stats.py
new file mode 100644
index 0000000..a919b7c
--- /dev/null
+++ b/mamotif/stats.py
@@ -0,0 +1,57 @@
+"""
+mamotif.stats
+-------------
+
+Statistical functions used in MAmotif.
+"""
+
+import numpy as np
+from scipy import stats
+
+
+def mamotif_t_test(m_values_pos, m_values_neg):
+    try:
+        t_stat, p_value = stats.ttest_ind(m_values_pos, m_values_neg,
+                                          equal_var=False)
+        if t_stat < 0:
+            p_right = 1 - p_value / 2
+        else:
+            p_right = p_value / 2
+        return t_stat, p_right
+    except:
+        return np.nan, np.nan
+
+
+def mamotif_ranksum_test(m_values_pos, m_values_neg):
+    try:
+        z_stat, p_value = stats.ranksums(m_values_pos, m_values_neg)
+        if z_stat < 0:
+            p_right = 1 - p_value / 2
+        else:
+            p_right = p_value / 2
+        return z_stat, p_right
+    except:
+        return np.nan, np.nan
+
+
+def adjust_p_values(p_values, correction='benjamin'):
+    n = len(p_values)
+    adjusted_p_values = []
+    if correction == 'benjamin':
+        order = np.argsort(p_values)
+        ranks = np.empty_like(order)
+        ranks[order] = np.arange(1, n + 1)
+        for p_value, rank in zip(p_values, ranks):
+            if np.isnan(p_value):
+                adjusted_p_values.append(np.nan)
+            else:
+                adjusted_p_values.append(min(1, p_value * n / rank))
+    elif correction == 'bonferroni':
+        for p_value in p_values:
+            if np.isnan(p_value):
+                adjusted_p_values.append(np.nan)
+            else:
+                adjusted_p_values.append(min(1, p_value * n))
+    else:
+        raise ValueError(f"invalid correction type: {correction}")
+    return adjusted_p_values