Repository: owkin/PyDESeq2 (python)

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)



Pull request:
BUG `beta` not initialized in `irls_solver` when the design is not a full rank matrix
#### Reference Issue or PRs
Fixes #132

#### What does your PR implement? Be specific.

* Initializes `beta` with log mean (normalized) counts when the design is not a full rank matrix.

This as such is not enough, because (as currently defined) the model cannot be fitted if the design matrix does not have full column rank. Hence, similarly to the [DESeq2](https://github.com/mikelove/DESeq2/blob/2c3c14cbfa09300e503e08302c60ddf112beacfd/R/AllClasses.R#L295) R package, this PR:

* Checks that the design matrix has full column rank at `DeseqDataSet` initialisation,
* Implements a test to check that a `ValueError` is thrown when the design matrix has deficient column rank.


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