Repository: owkin/PyDESeq2

Issue:
[BUG]` beta` referenced before assignment
**Describe the bug**
```
lib/python3.10/site-packages/pydeseq2/utils.py", line 471, in irls_solver
    mu = np.maximum(size_factors * np.exp(X @ beta), min_mu)
UnboundLocalError: local variable 'beta' referenced before assignment
```



![image](https://github.com/owkin/PyDESeq2/assets/39041821/c758f0ef-df80-4515-bfbd-b273f38b4fb6)



Discussion:

Patch:
diff --git a/pydeseq2/dds.py b/pydeseq2/dds.py
index 887f135..f611df6 100644
--- a/pydeseq2/dds.py
+++ b/pydeseq2/dds.py
@@ -217,6 +217,9 @@ class DeseqDataSet(ad.AnnData):
             intercept=True,
         )
 
+        # Check that the design matrix has full rank
+        self._check_full_rank_design()
+
         self.min_mu = min_mu
         self.min_disp = min_disp
         self.max_disp = np.maximum(max_disp, self.n_obs)
@@ -994,3 +997,15 @@ class DeseqDataSet(ad.AnnData):
 
         # Store normalized counts
         self.layers["normed_counts"] = self.X / self.obsm["size_factors"][:, None]
+
+    def _check_full_rank_design(self):
+        """Check that the design matrix has full column rank."""
+        rank = np.linalg.matrix_rank(self.obsm["design_matrix"])
+        num_vars = self.obsm["design_matrix"].shape[1]
+
+        if rank < num_vars:
+            raise ValueError(
+                "The design matrix is not full rank, so the model cannot be "
+                "fitted. Please remove the design variables that are linear "
+                "combinations of others."
+            )
diff --git a/pydeseq2/utils.py b/pydeseq2/utils.py
index 8b386fe..a58ea61 100644
--- a/pydeseq2/utils.py
+++ b/pydeseq2/utils.py
@@ -463,6 +463,9 @@ def irls_solver(
         y = np.log(counts / size_factors + 0.1)
         beta_init = solve(R, Q.T @ y)
         beta = beta_init
+    else:  # Initialise intercept with log base mean
+        beta = np.zeros(num_vars)
+        beta[0] = np.log(counts / size_factors).mean()
 
     dev = 1000.0
     dev_ratio = 1.0


Test patch:
diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py
index 37175e5..e60d0e3 100644
--- a/tests/test_edge_cases.py
+++ b/tests/test_edge_cases.py
@@ -116,7 +116,7 @@ def test_nan_factors():
 
 
 def test_one_factor():
-    """Test that a ValueError is thrown when the design factor takes only one value ."""
+    """Test that a ValueError is thrown when the design factor takes only one value."""
     counts_df = pd.DataFrame(
         {"gene1": [0, 1], "gene2": [4, 12]}, index=["sample1", "sample2"]
     )
@@ -126,6 +126,22 @@ def test_one_factor():
         DeseqDataSet(counts=counts_df, clinical=clinical_df, design_factors="condition")
 
 
+def test_rank_deficient_design():
+    """Test that a ValueError 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"]
+    )
+    clinical_df = pd.DataFrame(
+        {"condition": [0, 1], "batch": ["A", "B"]}, index=["sample1", "sample2"]
+    )
+
+    with pytest.raises(ValueError):
+        DeseqDataSet(
+            counts=counts_df, clinical=clinical_df, design_factors=["condition", "batch"]
+        )
+
+
 def test_reference_level():
     """Test that a ValueError is thrown when the reference level is not one of the
     design factor values."""