Repository: owkin/PyDESeq2 (python)

Issue:
ValueError: The first guess on the deviance function returned a nan.  This could be a boundary  problem and should be reported.
Hello, thank you for the great tool!

**Describe the bug**
When trying to run the deseq2 method, the following error arises: 
`
ValueError: The first guess on the deviance function returned a nan.  This could be a boundary  problem and should be reported.
`
**To Reproduce**
Version: `pydeseq2    0.4.1`
Code: 

```
#initialize DeseqDataSet
dds = DeseqDataSet(
    counts = counts, 
    metadata=pb.obs,
    design_factors = 'sample')

sc.pp.filter_genes(dds, min_cells = 1) #filter genes that are found in at least one sample

dds.deseq2() #call DESeq2
```
The last line produced the following output in addition to the ValueError:

> Fitting size factors...
> ... done in 0.00 seconds.
> 
> /home/jupyter/.local/lib/python3.8/site-packages/pydeseq2/utils.py:1014: RuntimeWarning: divide by zero encountered in divide
>   ((normed_counts - y_hat) ** 2 - y_hat) / ((num_samples - num_vars) * y_hat**2)
> /home/jupyter/.local/lib/python3.8/site-packages/pydeseq2/utils.py:1014: RuntimeWarning: invalid value encountered in divide
>   ((normed_counts - y_hat) ** 2 - y_hat) / ((num_samples - num_vars) * y_hat**2)
> Fitting dispersions...
> ... done in 2.34 seconds.
> 
> Fitting dispersion trend curve...
> /opt/conda/envs/NET_py_env/lib/python3.8/site-packages/statsmodels/genmod/families/family.py:177: RuntimeWarning: invalid value encountered in divide
>   return np.sum(resid_dev * freq_weights * var_weights / scale)

How should I proceed?

Pull request description:
BUG division by zero when ``num_samples == num_vars``
#### Reference Issue or PRs

Related to #174

#### What does your PR implement? Be specific.

Throws an error when the design matrix has the same number of rows and columns, as it is not possible to fit rough dispersions in this case. This reproduces the behavior of DESeq2.


Patch:
diff --git a/pydeseq2/utils.py b/pydeseq2/utils.py
index f456316..5b18559 100644
--- a/pydeseq2/utils.py
+++ b/pydeseq2/utils.py
@@ -1003,6 +1003,15 @@ def fit_rough_dispersions(
     """
 
     num_samples, num_vars = design_matrix.shape
+    # This method is only possible when num_samples > num_vars.
+    # If this is not the case, throw an error.
+    if num_samples == num_vars:
+        raise ValueError(
+            "The number of samples and the number of design variables are "
+            "equal, i.e., there are no replicates to estimate the "
+            "dispersion. Please use a design with fewer variables."
+        )
+
     # Exclude genes with all zeroes
     normed_counts = normed_counts[:, ~(normed_counts == 0).all(axis=0)]
 


Test patch:
diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py
index 1d6a036..5a2b8a7 100644
--- a/tests/test_edge_cases.py
+++ b/tests/test_edge_cases.py
@@ -163,8 +163,8 @@ def test_one_factor():
 
 
 def test_rank_deficient_design():
-    """Test that a ValueError is thrown when the design matrix does not have full column
-    rank."""
+    """Test that a UserWarning is thrown when the design matrix does not have full
+    column rank."""
     counts_df = pd.DataFrame(
         {"gene1": [0, 1], "gene2": [4, 12]}, index=["sample1", "sample2"]
     )
@@ -178,6 +178,28 @@ def test_rank_deficient_design():
         )
 
 
+def test_equal_num_vars_num_samples_design():
+    """Test that a ValueError is thrown when fitting dispersions if the design matrix
+    has eaual numbers of rows and columns."""
+    counts_df = pd.DataFrame(
+        {"gene1": [0, 1, 55], "gene2": [4, 12, 60]},
+        index=["sample1", "sample2", "sample3"],
+    )
+    metadata = pd.DataFrame(
+        {"condition": [0, 1, 0], "batch": ["A", "B", "B"]},
+        index=["sample1", "sample2", "sample3"],
+    )
+
+    dds = DeseqDataSet(
+        counts=counts_df, metadata=metadata, design_factors=["condition", "batch"]
+    )
+
+    dds.fit_size_factors()
+
+    with pytest.raises(ValueError):
+        dds.fit_genewise_dispersions()
+
+
 def test_reference_level():
     """Test that a ValueError is thrown when the reference level is not one of the
     design factor values."""