Moved stuff to libhts.

libhts/__init__.py

-from .libhts import do_deseq2
+from .libhts import do_deseq2, median_ratio_to_pseudo_ref_size_factors, set_de_status

libhts/libhts.py

 import pandas as pd
+# To compute geometric mean
+from scipy.stats.mstats import gmean
 from rpy2.robjects import r, pandas2ri, Formula, StrVector
 as_df = r("as.data.frame")
 from rpy2.robjects.packages import importr
@@ -38,3 +40,29 @@ def do_deseq2(cond_names, conditions, counts_data,
     res.index = counts_data.index
     return res
 
+
+def median_ratio_to_pseudo_ref_size_factors(counts_data):
+    """Adapted from DESeq paper (doi:10.1186/gb-2010-11-10-r106)"""
+    # Add pseudo-count to compute the geometric mean, then remove it
+    #pseudo_ref = (counts_data + 1).apply(gmean, axis=1) - 1
+    # Ignore lines with zeroes instead:
+    pseudo_ref = (counts_data[counts_data.prod(axis=1) > 0]).apply(gmean, axis=1)
+    def median_ratio_to_pseudo_ref(col):
+        return (col / pseudo_ref).median()
+    #size_factors = counts_data.apply(median_ratio_to_pseudo_ref, axis=0)
+    return counts_data[counts_data.prod(axis=1) > 0].apply(
+        median_ratio_to_pseudo_ref, axis=0)
+
+
+def set_de_status(row):
+    """Determines the up- or down-regulation status corresponding to a given
+    row of a deseq2 results table."""
+    if row["padj"] < 0.05:
+        #if row["log2FoldChange"] > 0:
+        if row["lfcMLE"] > 0:
+            return "up"
+        else:
+            return "down"
+    else:
+        return "NS"
+