Repository: owkin/PyDESeq2 (python)

Issue:
[BUG] `pydeseq2.ds.DeseqStats.lfc_shrink` gives wrong log2FoldChange signs
**Describe the bug**
When performing LFC shrinkage, the new version seems to flip the reference it uses. I've observed that the obtained `log2FoldChange` have opposite signs as the original contrast, which creates an inconsistency with the sign of `stat`.

**To Reproduce**
You can reproduce this bug with `pydeseq2 >= 0.3.1`. Here's a reproducible example:
```python
import pandas as pd
import numpy as np
from anndata import AnnData
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# Build toy data, note that G1 should be up in healthy, G2 does not change and G3 should be up in disease
var = pd.DataFrame(index=['G1', 'G2', 'G3'])
X = pd.DataFrame(
    [[16., 7., 12.],
     [22., 6., 1.],
     [10., 8., 30.],
     [3., 6., 20.],
     [1., 5., 17.]], columns=var.index
)
obs = pd.DataFrame(
    [['healthy'],
     ['healthy'],
     ['disease'],
     ['disease'],
     ['disease']], columns=['condition'], index=['S1', 'S2', 'S3', 'S4', 'S5']
    )
adata = AnnData(X, obs=obs, var=var, dtype=np.float32)

# Build DESeq2 object
dds = DeseqDataSet(
    adata=adata,
    design_factors='condition',
    refit_cooks=True,
    n_cpus=8,
)

# Compute LFCs
dds.deseq2()

# Extract contrast between disease vs healthy
stat_res = DeseqStats(dds, contrast=["condition", 'disease', 'healthy'], n_cpus=8)

# Compute Wald test
stat_res.summary()

# Store results before shrink
before_shrink = stat_res.results_df

# Assert that logFCs should have the same sign than stat
assert np.all(np.sign(before_shrink['log2FoldChange']) == np.sign(before_shrink['stat']))

# Shrink LFCs
stat_res.lfc_shrink(coeff='condition_healthy_vs_disease')

after_shrink = stat_res.results_df

# Assert that logFCs should have the same sign than stat
assert np.all(np.sign(after_shrink['log2FoldChange']) == np.sign(after_shrink['stat']))
```

**Expected behavior**
After LFC shrinkage, `log2FoldChange` should have the same sign than `stat`.

**Screenshots**
Before shrinkage:
![image](https://user-images.githubusercontent.com/44896790/234716907-9d16716b-b864-47f3-a22f-06e867ccaefc.png)

After shrinkage:
![image](https://user-images.githubusercontent.com/44896790/234716974-68b6b232-0926-47c5-a221-dc76ea73f8f5.png)

**Desktop (please complete the following information):**
 - OS: linux

**Additional context**
May I ask why now the method `pydeseq2.ds.DeseqStats.lfc_shrink` requires `coeff`? Before (between version `0.3.0` and `0.3.1`) it could be run based on the contrast specified already in the `DeseqStats` object. Now you need to manually check `stat_res.LFC` for the name which feels weird.

I'm aware that if you provide `ref_level=['condition', 'healthy']` to the previous `DeseqDataSet` it will work but then there are two times where the contrast needs to be specified and if there are discrepancies, like here, errors will rise. IMO it makes more sense to specify the contrast at the `DeseqStats` level instead of the `DeseqStats` one.


Pull request description:
BUG wrong `lfc_shrink` sign display
#### Reference Issue or PRs
Fixes #121 


#### What does your PR implement? Be specific.
This PR addresses #121 (wrong display of shrunk factors). More precisely:
- the message output by `lfc_shrink` was fixed to match the coefficient that was actually shrunk
- when the coefficient is not available, the error message now explains how to set `ref_level` in `DeseqDataSet`
- it is now possible to set `coeff=None`
- `coeff=None` is tested in `test_edge_cases.py`


Patch:
diff --git a/README.md b/README.md
index c521e40..825f61b 100644
--- a/README.md
+++ b/README.md
@@ -105,7 +105,7 @@ documentation to see how you can contribute to PyDESeq2.
 
 ### 2 - Create a conda environment
 
-Run `conda env create -n pydeseq2 python=3.8` (or higher python version) to create the `pydeseq2` environment and then activate it:
+Run `conda create -n pydeseq2 python=3.8` (or higher python version) to create the `pydeseq2` environment and then activate it:
 `conda activate pydeseq2`.
 
 `cd` to the root of the repo and run `pip install -e ."[dev]"` to install in developer mode.
diff --git a/pydeseq2/dds.py b/pydeseq2/dds.py
index 7ae1a35..f611df6 100644
--- a/pydeseq2/dds.py
+++ b/pydeseq2/dds.py
@@ -16,6 +16,7 @@ from joblib import parallel_backend
 from scipy.optimize import minimize
 from scipy.special import polygamma  # type: ignore
 from scipy.stats import f  # type: ignore
+from scipy.stats import trim_mean  # type: ignore
 from statsmodels.tools.sm_exceptions import DomainWarning  # type: ignore
 
 from pydeseq2.preprocessing import deseq2_norm
@@ -216,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)
@@ -226,7 +230,11 @@ class DeseqDataSet(ad.AnnData):
         self.batch_size = batch_size
         self.joblib_verbosity = joblib_verbosity
 
-    def vst(self, use_design: bool = False) -> None:
+    def vst(
+        self,
+        use_design: bool = False,
+        fit_type: Literal["parametric", "mean"] = "parametric",
+    ) -> None:
         """
         Fit a variance stabilizing transformation, and apply it to normalized counts.
 
@@ -237,6 +245,11 @@ class DeseqDataSet(ad.AnnData):
         use_design : bool
             Whether to use the full design matrix to fit dispersions and the trend curve.
             If False, only an intercept is used. (default: ``False``).
+        fit_type : str
+            Either "parametric" or "mean" for the type of fitting of dispersions to the
+            mean intensity. parametric - fit a dispersion-mean relation via a robust
+            gamma-family GLM. mean - use the mean of gene-wise dispersion estimates.
+            (default: ``"parametric"``).
         """
 
         # Start by fitting median-of-ratio size factors, if not already present.
@@ -256,19 +269,34 @@ class DeseqDataSet(ad.AnnData):
             )
             # Fit the trend curve with an intercept design
             self.fit_genewise_dispersions()
-            self.fit_dispersion_trend()
+            if fit_type == "parametric":
+                self.fit_dispersion_trend()
 
             # Restore the design matrix and free buffer
             self.obsm["design_matrix"] = self.obsm["design_matrix_buffer"].copy()
             del self.obsm["design_matrix_buffer"]
 
         # Apply VST
-        a0, a1 = self.uns["trend_coeffs"]
-        cts = self.layers["normed_counts"]
-        self.layers["vst_counts"] = np.log2(
-            (1 + a1 + 2 * a0 * cts + 2 * np.sqrt(a0 * cts * (1 + a1 + a0 * cts)))
-            / (4 * a0)
-        )
+        if fit_type == "parametric":
+            a0, a1 = self.uns["trend_coeffs"]
+            cts = self.layers["normed_counts"]
+            self.layers["vst_counts"] = np.log2(
+                (1 + a1 + 2 * a0 * cts + 2 * np.sqrt(a0 * cts * (1 + a1 + a0 * cts)))
+                / (4 * a0)
+            )
+        elif fit_type == "mean":
+            gene_dispersions = self.varm["genewise_dispersions"]
+            use_for_mean = gene_dispersions > 10 * self.min_disp
+            mean_disp = trim_mean(gene_dispersions[use_for_mean], proportiontocut=0.001)
+            self.layers["vst_counts"] = (
+                2 * np.arcsinh(np.sqrt(mean_disp * self.layers["normed_counts"]))
+                - np.log(mean_disp)
+                - np.log(4)
+            ) / np.log(2)
+        else:
+            raise NotImplementedError(
+                f"Found fit_type '{fit_type}'. Expected 'parametric' or 'mean'."
+            )
 
     def deseq2(self) -> None:
         """Perform dispersion and log fold-change (LFC) estimation.
@@ -969,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/ds.py b/pydeseq2/ds.py
index bb8b105..3cff5f1 100644
--- a/pydeseq2/ds.py
+++ b/pydeseq2/ds.py
@@ -267,25 +267,52 @@ class DeseqStats:
             self.statistics.loc[self.dds.new_all_zeroes_genes] = 0.0
             self.p_values.loc[self.dds.new_all_zeroes_genes] = 1.0
 
-    def lfc_shrink(self, coeff: str) -> None:
+    def lfc_shrink(self, coeff: Optional[str] = None) -> None:
         """LFC shrinkage with an apeGLM prior :cite:p:`DeseqStats-zhu2019heavy`.
 
         Shrinks LFCs using a heavy-tailed Cauchy prior, leaving p-values unchanged.
 
         Parameters
         ----------
-        coeff : str
-            The LFC coefficient to shrink.
+        coeff : str or None
+            The LFC coefficient to shrink. If set to ``None``, the method will try to
+            shrink the coefficient corresponding to the ``contrast`` attribute.
             If the desired coefficient is not available, it may be set from the
             :class:`pydeseq2.dds.DeseqDataSet` argument ``ref_level``.
+            (default: ``None``).
         """
 
-        if coeff not in self.LFC.columns:
+        contrast_level = f"{self.contrast[0]}_{self.contrast[1]}_vs_{self.contrast[2]}"
+
+        if coeff is not None:
+            if coeff not in self.LFC.columns:
+                split_coeff = coeff.split("_")
+                if len(split_coeff) == 4:
+                    raise KeyError(
+                        f"The coeff argument '{coeff}' should be one the LFC columns. "
+                        f"The available LFC coeffs are {self.LFC.columns[1:]}. "
+                        f"If the desired coefficient is not available, please set "
+                        f"`ref_level = [{split_coeff[0]}, {split_coeff[3]}]` "
+                        f"in DeseqDataSet and rerun."
+                    )
+                else:
+                    raise KeyError(
+                        f"The coeff argument '{coeff}' should be one the LFC columns. "
+                        f"The available LFC coeffs are {self.LFC.columns[1:]}. "
+                        f"If the desired coefficient is not available, please set the "
+                        f"appropriate`ref_level` in DeseqDataSet and rerun."
+                    )
+        elif contrast_level not in self.LFC.columns:
             raise KeyError(
-                f"The coeff argument ('{coeff}') should be one the LFC columns. "
-                f"If not available,  it can be set from DeseqDataSet's `ref_level` "
-                f"argument."
+                f"lfc_shrink's coeff argument was set to None, but the coefficient "
+                f"corresponding to the contrast {self.contrast} is not available."
+                f"The available LFC coeffs are {self.LFC.columns[1:]}. "
+                f"If the desired coefficient is not available, please set "
+                f"`ref_level = [{self.contrast[0]}, {self.contrast[2]}]` "
+                f"in DeseqDataSet and rerun."
             )
+        else:
+            coeff = contrast_level
 
         coeff_idx = self.LFC.columns.get_loc(coeff)
 
@@ -356,9 +383,14 @@ class DeseqStats:
             self.results_df["log2FoldChange"] = self.LFC.iloc[:, coeff_idx] / np.log(2)
             self.results_df["lfcSE"] = self.SE / np.log(2)
 
+            # Get the corrresponding factor, tested and reference levels of the shrunk
+            # coefficient
+            split_coeff = coeff.split("_")
+            # coeffs are of the form "factor_A_vs_B", hence "factor" is split_coeff[0],
+            # "A" is split_coeff[1] and "B" split_coeff[3]
             print(
                 f"Shrunk Log2 fold change & Wald test p-value: "
-                f"{self.contrast[0]} {self.contrast[1]} vs {self.contrast[2]}"
+                f"{split_coeff[0]} {split_coeff[1]} vs {split_coeff[3]}"
             )
 
             display(self.results_df)
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/data/single_factor/r_mean_vst.csv b/tests/data/single_factor/r_mean_vst.csv
new file mode 100644
index 0000000..b75ee76
--- /dev/null
+++ b/tests/data/single_factor/r_mean_vst.csv
@@ -0,0 +1,11 @@
+"","sample1","sample2","sample3","sample4","sample5","sample6","sample7","sample8","sample9","sample10","sample11","sample12","sample13","sample14","sample15","sample16","sample17","sample18","sample19","sample20","sample21","sample22","sample23","sample24","sample25","sample26","sample27","sample28","sample29","sample30","sample31","sample32","sample33","sample34","sample35","sample36","sample37","sample38","sample39","sample40","sample41","sample42","sample43","sample44","sample45","sample46","sample47","sample48","sample49","sample50","sample51","sample52","sample53","sample54","sample55","sample56","sample57","sample58","sample59","sample60","sample61","sample62","sample63","sample64","sample65","sample66","sample67","sample68","sample69","sample70","sample71","sample72","sample73","sample74","sample75","sample76","sample77","sample78","sample79","sample80","sample81","sample82","sample83","sample84","sample85","sample86","sample87","sample88","sample89","sample90","sample91","sample92","sample93","sample94","sample95","sample96","sample97","sample98","sample99","sample100"
+"gene1",3.43178816207728,0.794238482935036,2.3243923299861,0.929900166549763,1.08424327046672,-0.92587472327995,-0.92587472327995,1.42407402454225,3.06261011677541,3.89776446507256,1.94554812217562,2.53649012674806,1.84546698268671,4.88837163614702,2.96427801123802,5.42291256184733,0.710106671447233,5.02390639115651,2.04078009680906,0.727745958147907,1.69643682198544,3.26381753389741,1.46416434886128,4.04092013780374,0.792459761890463,4.56257204314427,-0.92587472327995,3.34060543576917,1.02517707019824,0.743086416612365,2.85232029689219,-0.92587472327995,3.23209472527434,2.34222356686588,1.73406037037047,1.51532010962472,1.93351221754072,-0.92587472327995,1.56254092259907,3.74420129919689,1.11817213179645,3.4452844897408,1.69123555747773,4.4894136364452,0.84929877436974,-0.92587472327995,3.59551957296111,1.52323842567228,4.27420396999495,2.11767998712677,3.76683338455236,3.00839774523434,4.00835989325415,3.62309553308489,2.3533340818367,1.83004272995896,2.80603404546218,4.72284250804227,3.84962648983205,1.93275360917771,2.59097031968391,0.999863563155015,3.12565473902118,4.23943233009275,2.69889333519741,3.70272392373991,4.17508121870462,3.0062573294081,1.50588316223104,4.79480558779146,4.8419088744785,3.11538829035537,2.64689547542903,3.08588631633697,3.24695150951301,2.78165423016523,2.20656139405463,3.43963012628544,3.26883038578985,4.39094508565612,1.9247951108824,1.64287408542231,3.09478616745425,3.51600646472837,4.83158419828075,2.6557879676613,3.80015136860503,3.27614123410574,3.72271204291999,2.84032013497841,-0.92587472327995,3.21969683221821,2.8892267539958,4.82543539609511,1.42425597620222,2.70594576764465,1.18588608349968,3.86288947482366,4.30820361384928,5.090795405629
+"gene2",4.17989032061667,5.25014198336425,2.3243923299861,3.46022075911297,3.82919533080377,3.04332194900246,4.05252963322757,4.13005732295277,4.85385369224168,3.24799621688382,4.26465682505567,4.33815210000116,3.14986135219646,3.99489829450734,5.38435518383919,4.07663842163417,4.4667493519774,3.84361366379629,3.69182304501022,3.89776446507256,2.97060519543833,4.39661725186723,4.44740951730844,4.32859490529157,4.01277412047786,4.44274003513364,4.54627463710865,4.61200266534663,4.15513753984006,5.02266095230433,4.34209536362991,4.63804993367197,3.94378658405492,3.19660760402319,2.70271837156868,3.74420129919689,3.33802483932071,3.70639162532982,4.29696200256275,4.63382673715581,4.83019510457191,4.27499356874843,4.5646855456792,3.83737109230778,2.46359963432519,4.20800929731357,4.60073348753461,3.85396714776233,4.98011588766161,4.51238759537133,4.96511752850299,4.50248209944965,5.36454683623403,4.56214118498245,4.49991152521692,5.51659544078661,4.71661404870768,4.86919024759426,4.40011133967893,4.878729719055,4.00812849355124,4.92389041130604,4.12219881004377,3.1425599213058,3.40896458785202,4.35907120179453,4.41002398133365,3.73452980721184,4.52545703574313,3.5950473352912,2.98671773279319,4.51630255121467,5.20759220595933,5.54280248254242,4.8889742741107,4.02520898664886,3.3848930362253,2.98068877989954,5.29767153666374,4.89519384977852,4.33593536476715,4.80883307032224,4.24070413553716,4.45041313844967,2.90111578094191,6.56990345084548,5.41293821038467,5.34885504360545,5.86417085093878,3.86381434376906,5.04648271670275,5.34581074405528,4.14240440341942,3.71913262157252,4.4496479666486,4.42778225057847,4.31150477843712,5.63528431845285,3.96674161996297,3.92523323836747
+"gene3",2.08479263722626,1.39493931723213,3.59220466704363,1.56298564962319,3.03707223492813,1.46387004136244,4.31685581756083,3.40097533517357,3.83154710207773,2.84381169647346,3.75203791032984,2.53649012674806,2.16403191147716,3.61012894658851,3.86868743641941,3.16028205412776,3.21740061691931,2.20987503397319,-0.92587472327995,2.27760254323782,0.715366740834192,1.24235125406151,2.46449932413447,0.85335032838488,1.80370356295615,0.882583251221189,4.41510745579085,3.34060543576917,3.31815624558574,-0.92587472327995,1.1152182414957,-0.92587472327995,2.80609239452915,2.01388555561562,2.86801727649087,1.93911444315392,1.93351221754072,3.80528596384481,0.620544709358,3.74420129919689,3.8812569824696,3.4452844897408,1.99999829583264,3.028445311188,2.68385546304285,-0.92587472327995,3.72201679699642,1.14130215461576,3.06071264648486,1.39100700150718,2.25356147506667,-0.92587472327995,5.27386098038003,1.12447390451642,0.954284128957104,-0.92587472327995,1.56187672318477,-0.92587472327995,1.520016035625,2.25650097635248,0.93465676192635,1.36373732131704,2.5415239514582,2.11628869905417,2.69889333519741,2.2451590724531,2.06315755917287,0.826661027029581,1.50588316223104,0.865401542930877,2.36295062791096,3.92913742084781,2.83681672993213,0.704686068691239,2.05742424858567,1.06567840154477,0.682281023814083,2.30030550388748,-0.92587472327995,0.753619435479918,2.30467238424358,2.6795702808216,2.51248482725021,0.88706342875241,1.35502464654549,1.85711327033525,3.26636983696285,2.2023583090569,2.21543652487421,1.10676819682919,3.63068792109372,2.22628011499165,4.03386189839199,-0.92587472327995,0.817770391578745,1.73681897398091,2.33053242601713,1.87993907608618,2.04055190161778,2.51688846346766
+"gene4",6.7391571542759,5.75604606807718,7.50424360070479,6.6558479022771,7.34704132668137,5.06849318701028,7.47360382634241,5.72971342245206,6.61287455142089,6.93119637116649,7.12740623137526,5.16964691115809,6.81703962234509,5.92316954927647,7.45330372957535,6.77668566065011,7.25802924059534,7.37304297109047,7.39364954135367,5.69790553156317,6.39952225388508,6.46422961158066,6.75016827604946,6.03272129211965,5.96149593633252,7.57165780503787,7.32520517759504,7.28387317965898,7.0064443254571,6.79660615303825,7.00336851962609,5.94627292216031,6.46546528635387,6.3840947526446,7.99790261112106,6.72390211209415,7.44046842757332,6.72179844637657,6.6166334323622,6.72390211209415,7.35955968020412,6.63957479424225,7.07987475624679,7.36879717840363,7.04373563295402,5.8515028143795,6.62004696801525,6.84364391388558,6.2845626693185,6.85579664584554,7.01530733390636,5.81233875166987,6.66438568992049,5.78121292352769,6.89793005324003,5.21625119182065,6.84146449481982,5.20861471199194,6.4283717036607,6.4292788422282,6.1959864345532,3.92101882019935,7.02396010864191,6.64879300601062,6.59994380769147,6.08190000609387,6.49891675413818,5.96446676363181,6.30295891874055,6.35957062861121,7.06287532844036,4.98301554636892,6.71792563847355,6.2574332632185,6.55796419853726,5.03164665746839,6.76840626059335,5.59086324689928,6.62472685912832,5.63523226742711,6.85079506716932,6.67971755785652,6.68102807005755,6.89920066383141,6.98767837847414,6.70413072341641,6.85417285428167,6.04493387555856,6.92526996951615,6.75516865141753,6.95306167035072,6.32719889726835,7.17255845844929,6.50893687720827,6.00785383889893,5.75034163058625,6.59216111798675,6.62004696801525,6.12678848319447,6.06859918902884
+"gene5",3.97097653081952,3.36162822476565,4.46588827648748,5.49030773049257,4.33702544608528,4.61599777273581,3.72850832640206,4.13005732295277,3.49794087954113,4.80811537613037,2.27003371474756,4.79533649961694,5.02063415132218,4.54863061274862,5.03199083296529,4.75807910804578,3.86589904014313,2.82158604787481,5.04277923591031,5.18958552443388,3.10428582507447,4.18587989873431,4.90454061682071,4.39235525955955,5.30931616741387,2.25078283907753,3.9298654122258,6.25719824820746,4.06149829139974,3.92529388514136,3.36971724605468,4.81377690367881,4.3346719885916,4.83461406417124,4.75336568179647,5.63136706523312,4.38660508777375,5.24933006417711,4.29696200256275,4.1052246069752,4.20762178158816,4.16586014070359,4.51793840610849,4.27797072707265,3.32941033817102,4.41892014658786,4.60073348753461,4.80193956928288,3.72672336340452,3.23280318107341,4.22108024946275,4.5365075264258,4.56237746641135,4.40004587180739,5.24979589212791,5.48922397213177,3.84029715352425,5.4602653789638,5.36437309904432,5.4149984561965,4.96132819780113,4.92389041130604,5.07892817028401,5.10911750011833,5.25003620057548,5.93059434499213,4.33588679439673,5.37124845212033,5.2522314869794,5.3123778841454,4.02105956465063,4.44571043101009,3.15439623070099,3.77953385451771,3.98712207357607,5.61858088121626,3.3848930362253,5.56115258820145,4.03128416354184,5.23786698058637,4.62672854618806,5.78186894365141,5.16344644761721,4.51158970232,5.99952305005083,4.38465475594441,6.18391001039403,4.54976616167457,4.17601274515359,5.37924774938093,5.71879893410327,5.27889478679021,4.80193956928288,6.46004642101066,4.80343485127579,3.27648483902399,5.50467044370666,3.9911917915754,5.81559154179159,4.42555810888592
+"gene6",-0.92587472327995,3.2380819437829,1.99699257177663,3.32274368244967,1.75140184629963,1.88239230210189,4.46947388881781,1.42407402454225,2.43365400589861,-0.92587472327995,3.11718246368632,-0.92587472327995,1.43044973338473,3.28105370334591,2.96427801123802,-0.92587472327995,3.21740061691931,1.0936069856165,3.62007445393761,1.31171582666095,3.33919100713179,4.06769434185303,3.04370306286401,-0.92587472327995,3.67589846072944,1.50463705270734,2.60426544539576,4.75790693549745,-0.92587472327995,1.33096705436949,3.28011806464201,1.14886719583114,2.19556815410358,1.58334685584163,2.04562964523027,2.26322949706866,1.40710392455626,3.60019047088059,1.56254092259907,1.51532010962472,2.24215952798045,1.228884131548,1.69123555747773,3.83737109230778,1.46342406364546,-0.92587472327995,4.04593320540057,2.8892267539958,3.06071264648486,0.791086417446387,3.47140163324332,0.996498578864634,3.29395937315881,-0.92587472327995,2.84427388789767,2.98175326494718,2.80603404546218,2.87896552930181,2.26869706505936,-0.92587472327995,2.32541074277567,-0.92587472327995,2.76388198988043,3.73379291131391,4.60549567490039,-0.92587472327995,4.08740301404419,3.91335864664364,1.12585760011209,2.22598279786236,-0.92587472327995,2.09232366410416,1.84926961596114,2.45538161923126,3.07572827208226,1.72889046562596,0.682281023814083,1.97418698055839,3.64579508151224,1.74983713863031,0.507101031463376,1.64287408542231,1.92635726857939,3.39111437639062,2.0752524566744,1.17187235749091,-0.92587472327995,2.7628938947394,3.42806106663339,-0.92587472327995,3.50487669912145,-0.92587472327995,2.272447702158,-0.92587472327995,3.51527418442247,2.04856676848526,2.95184190310581,1.87993907608618,3.64411938926023,4.42555810888592
+"gene7",3.81269559865565,5.90526200245796,4.86399796209953,4.85712889803863,5.29970059406502,5.06849318701028,5.4719880805878,5.19341940322457,4.32988702315417,4.63987077564419,4.79904613937181,4.07014090620947,4.51461533881884,4.82652878410092,3.48654606575442,4.56466872061655,4.73494579868074,3.42172751107245,4.92879470839935,4.3975201256655,5.00534524569763,4.74342784372977,3.67522748721101,5.44924376944222,3.67589846072944,4.61894764930429,4.34475659320036,4.4496479666486,4.40427406340524,3.92529388514136,3.45404998264196,4.56128691782555,5.11007101658118,5.34891370643955,4.62136974471134,3.93600340561755,4.96088690931289,3.89781703245662,4.95261110670347,4.39363209151012,3.8812569824696,6.83032392105942,3.94141048990625,4.42232205685277,5.25003214726181,5.16721870619633,5.04934821228878,4.50957799651206,5.86843713619107,4.28313509341096,3.89452282897172,4.8387520669373,5.44986621198554,3.8908717069195,4.0403783500082,4.25572430360428,4.56278893890849,3.70145828336292,4.89339937187823,4.878729719055,5.36430601987849,4.99930108081322,3.95284945294688,4.6772932302417,3.88180041548788,4.41084284594692,5.70400303076413,3.82672093566504,4.84232807381617,3.3558119653633,6.07412426340524,6.00889589342914,4.90420880993987,6.02891940055087,4.66636838536957,4.86741695559438,5.5452608009932,4.78473465169317,3.94415580222581,5.53950452156322,5.12307422583551,4.19444512089503,4.82342748891026,4.31966960701174,4.35066408898929,5.75432771866845,4.18856791730183,4.40969615188884,5.03067081578386,5.6853489348002,2.35844585417078,4.28297499869708,5.34642612255901,4.26641949907168,4.61200266534663,5.05936293931835,4.73226921952812,6.49591126963135,5.10504712560234,4.95101438336282
+"gene8",5.49138515105038,4.80535998210284,5.13502353897528,4.08763066327762,5.16160983412869,5.28449871792694,4.67180349323477,6.18997720766278,5.23732320035939,5.36171805499149,5.33211648477166,5.24979589212791,4.67113041868971,6.35406333243083,3.33266821356054,6.19859142560841,3.32996903471086,5.7639211895824,5.0699210986472,5.60957460083053,5.58610243159403,5.2566807574007,4.32261719059448,6.87267128335333,5.58633543714352,6.57247002608383,6.13032375238654,4.89039195887469,6.40847110075396,5.05716282373857,5.55398968889593,5.94627292216031,5.38250492679738,4.19745740759325,6.03294575129939,5.54579720532638,5.0377489581034,5.93450271516517,4.64266553561399,5.83746519946583,5.70371479214432,5.57152770044192,5.38437233655922,5.58407403386564,5.34922409457279,5.38364185260414,5.23027343688977,5.43438784918505,4.52441920617362,5.47930916158157,4.84392395711782,5.48668065830373,4.21760973531583,5.71383363347951,5.08477982759049,4.49480458586405,3.84029715352425,4.86919024759426,5.1864159384297,4.24869509509713,5.62125391475239,4.99930108081322,5.47244220782269,5.01267282419977,5.25003620057548,5.99980738214571,5.67424288060765,5.94405769123258,4.31375141193251,5.69246295103871,5.20834511069004,4.76846428281604,3.91138731058602,5.58832181967136,5.29133873555656,5.1315835760718,5.98560598454242,5.90653694440739,5.22567429793846,4.6875205026063,2.95941545579286,4.96661645063311,4.56131882081547,5.47849307859205,5.80503995047424,5.14772825939178,4.66696340190874,5.44764081670182,3.96713396781713,5.09771589052546,4.9044759915309,4.932137498724,4.86662177538351,5.33963258498094,4.89039195887469,5.15818045523617,5.84548164340383,5.23027343688977,4.38226770842211,6.95126529517045
+"gene9",5.35460669625138,5.61458136669515,6.04940409777487,5.08633861318406,5.25512354757598,5.28449871792694,5.43540288965987,5.35497156370465,6.27306660930164,4.8854050912463,2.53355064556866,5.27555200421566,5.87394679929587,5.47666286380591,4.75807910804578,5.03199083296529,5.15634461965484,6.05423725325639,4.35337766916684,5.60957460083053,5.33834132053998,5.71521948139777,5.02954291061152,5.15068356833244,5.50866924935231,4.61894764930429,4.92935764737061,4.26665075304699,5.32903264333853,5.36319062105616,5.13687754220615,4.74603439453723,6.48486318959793,6.43591660271651,4.52613974751239,5.01973314834636,5.44761551019176,5.0217900841555,5.46817247894165,5.01973314834636,3.99846012945656,5.06433276875761,4.99835607406753,5.70185968999486,5.06785711197069,6.5997722332741,5.00033321553451,4.98788199593517,5.45192470626727,4.96363500118135,6.10810725825648,4.01557865837787,4.31183864253693,5.40759633943552,3.74056021588874,5.64616433735923,5.330961865648,5.54979575296132,4.58394398112198,5.01161175775046,4.09550877539301,4.06848484314728,2.76388198988043,5.61529874232569,6.3318583100297,5.11953259811092,4.54775184425308,5.76947309851899,5.04815984576871,5.37777595196735,5.20834511069004,5.37355541843204,5.02557718997444,4.06388842325048,5.29133873555656,5.47378302845993,4.97624702934358,5.1061145889738,5.22567429793846,4.1527659364032,4.16508718145026,5.0629623578111,4.44149164346818,4.92528256420026,4.10203679773782,4.71064588331885,4.89260308960655,5.24280200969325,4.85973916529084,4.81213729313319,4.04627303651183,5.4097591052487,4.98788199593517,6.20703539980641,5.88621264309139,5.44584315528753,5.11465705709758,4.50243744437511,6.38261881966225,5.59211194734463
+"gene10",1.77086401486503,3.10287657553229,1.56795736257554,3.32274368244967,2.80822497787611,4.0346051864541,2.94187130330695,3.51502727273954,2.43365400589861,3.12552582459743,1.94554812217562,3.09730081762681,2.64258250800084,4.76191377831068,1.69122289759055,3.86868743641941,2.96165999986297,3.84361366379629,3.0926374003611,2.84381169647346,2.25834106748151,5.02277138332688,3.95314279875382,3.87219734612093,3.23510380349769,2.7351891955462,3.0217871883997,4.95233224265644,4.32592003832783,1.33096705436949,-0.92587472327995,2.89989254743684,5.53917755992322,3.02609617741563,2.70271837156868,2.94044849412914,4.22061485853965,3.70639162532982,2.93846553358845,3.26100557906343,-0.92587472327995,3.7778484078723,1.99999829583264,3.83737109230778,3.45377539162844,3.81925078546313,3.94584867790651,2.8892267539958,3.25834086843592,3.35630068223929,1.80293619815936,3.86293298678138,1.80022735376905,5.36430601987849,3.20880196045896,3.5047504942015,3.62371732095961,4.93714575892283,4.84621966863621,1.93275360917771,4.46708730052684,2.08553831426531,2.95618399685314,3.95776024819014,4.01021894768448,3.86011287888135,2.06315755917287,4.57560492351408,1.80171484853465,5.27853206527517,2.98671773279319,4.44571043101009,4.77167689792595,2.80520653630052,2.05742424858567,3.20668599868345,1.25448031766566,3.15052944014657,4.89525265182405,2.06242371345113,1.39941250626637,3.10038236969977,3.92029146284768,-0.92587472327995,1.35502464654549,4.7817993595369,3.64327175634213,2.90947839793305,3.72271204291999,5.61469620434808,3.85343330896099,0.865607163639861,5.34642612255901,4.87602280495839,2.15676386849244,0.744270982466206,2.67484335644609,3.72201679699642,3.64411938926023,3.15166291979988
diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py
index 37175e5..4fe663b 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."""
@@ -145,7 +161,7 @@ def test_reference_level():
 
 def test_lfc_shrinkage_coeff():
     """Test that a KeyError is thrown when attempting to shrink an unexisting LFC
-    coefficient.
+    coefficient, but that setting coeff=None runs bug-free.
     """
 
     counts_df = load_example_data(
@@ -165,10 +181,22 @@ def test_lfc_shrinkage_coeff():
     )
     dds.deseq2()
 
-    res = DeseqStats(dds)
-    res.summary()
+    # Check that coeff=None gives the same results as the default
+    # coefficient condition_B_vs_A
+    res_1 = DeseqStats(dds)
+    res_1.summary()
+    res_1.lfc_shrink(coeff=None)
+    shrunk_lfcs_1 = res_1.results_df["log2FoldChange"].copy()
+
+    res_2 = DeseqStats(dds)
+    res_2.summary()
+    res_2.lfc_shrink(coeff="condition_B_vs_A")
+    shrunk_lfcs_2 = res_2.results_df["log2FoldChange"].copy()
+
+    assert (shrunk_lfcs_1 == shrunk_lfcs_2).all()
+
     with pytest.raises(KeyError):
-        res.lfc_shrink(coeff="this_coeff_does_not_exist")
+        res_1.lfc_shrink(coeff="this_coeff_does_not_exist")
 
 
 def test_indexes():
diff --git a/tests/test_pydeseq2.py b/tests/test_pydeseq2.py
index ba33114..0e6fcc2 100644
--- a/tests/test_pydeseq2.py
+++ b/tests/test_pydeseq2.py
@@ -449,6 +449,39 @@ def test_vst(tol=0.02):
     ).max().max() < tol
 
 
+def test_mean_vst(tol=0.02):
+    """
+    Test the output of VST with ``fitType="mean"`` compared with DESeq2.
+    """
+
+    # Load data
+    counts_df = load_example_data(
+        modality="raw_counts",
+        dataset="synthetic",
+        debug=False,
+    )
+
+    clinical_df = load_example_data(
+        modality="clinical",
+        dataset="synthetic",
+        debug=False,
+    )
+
+    # Load R data
+    test_path = str(Path(os.path.realpath(tests.__file__)).parent.resolve())
+
+    r_vst = pd.read_csv(
+        os.path.join(test_path, "data/single_factor/r_mean_vst.csv"), index_col=0
+    ).T
+
+    # Test blind design
+    dds = DeseqDataSet(
+        counts=counts_df, clinical=clinical_df, design_factors=["condition"]
+    )
+    dds.vst(use_design=False, fit_type="mean")
+    assert (np.abs(r_vst - dds.layers["vst_counts"]) / r_vst).max().max() < tol
+
+
 def test_ref_level():
     """Test that DeseqDataSet columns are created according to the passed reference
     level, if any.