More error-handling in do_deseq2.

libhts/libhts.py

+import warnings
+
+
+def formatwarning(message, category, filename, lineno, line):
+    """Used to format warning messages."""
+    return "%s:%s: %s: %s\n" % (filename, lineno, category.__name__, message)
+
+
+warnings.formatwarning = formatwarning
+
 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.rinterface import RRuntimeError
 from rpy2.robjects.packages import importr
 deseq2 = importr("DESeq2")
 
+
 def do_deseq2(cond_names, conditions, counts_data,
               formula=None, contrast=None, deseq2_args=None):
     """Runs a DESeq2 differential expression analysis."""
@@ -31,7 +43,39 @@ def do_deseq2(cond_names, conditions, counts_data,
         countData=counts_data,
         colData=col_data,
         design=formula)
-    dds = deseq2.DESeq(dds, betaPrior=deseq2_args["betaPrior"])
+    # dds = deseq2.DESeq(dds, betaPrior=deseq2_args["betaPrior"])
+    # Decompose into the 3 steps to have more control on the options
+    try:
+        dds = deseq2.estimateSizeFactors_DESeqDataSet(dds, type="ratio")
+    except RRuntimeError as e:
+        if sum(counts_data.prod(axis=1)) == 0:
+            msg = "".join(["Error occurred in estimateSizeFactors:\n%s\n" % e,
+                           "This is probably because every gene has at least one zero.\n",
+                           "We will try to use the \"poscounts\" method instead."])
+            warnings.warn(msg)
+            dds = deseq2.estimateSizeFactors_DESeqDataSet(dds, type="poscounts")
+        else:
+            raise
+    try:
+        dds = deseq2.estimateDispersions_DESeqDataSet(dds, fitType="parametric")
+    except RRuntimeError as e:
+        msg = "".join(["Error occurred in estimateDispersions:\n%s\n" % e,
+                    "We will try with fitType=\"local\"."])
+        warnings.warn(msg)
+        try:
+            dds = deseq2.estimateDispersions_DESeqDataSet(dds, fitType="local")
+        except RRuntimeError as e:
+            msg = "".join(["Error occurred in estimateDispersions:\n%s\n" % e,
+                        "We will try with fitType=\"mean\"."])
+            warnings.warn(msg)
+            try:
+                dds = deseq2.estimateDispersions_DESeqDataSet(dds, fitType="mean")
+            except RRuntimeError as e:
+                msg = "".join(["Error occurred in estimateDispersions:\n%s\n" % e,
+                            "We give up."])
+                warnings.warn(msg)
+                raise
+    dds = deseq2.nbinomWaldTest(dds, betaPrior=deseq2_args["betaPrior"])
     res = pandas2ri.ri2py(as_df(deseq2.results(
         dds,
         contrast=contrast,