Repository: owkin/PyDESeq2

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?

Discussion:

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."""