Repository: owkin/PyDESeq2

Issue:
Implement fitType="mean" for vst
DESeq2 has a fit type "mean" for the vst which is used as the fallback when the parametric fit fails. It would be nice to have an equivalent in the python implementation.

Here's my attempt based on https://github.com/mikelove/DESeq2/blob/devel/R/vst.R#L188 and https://github.com/mikelove/DESeq2/blob/devel/R/core.R#L892

```
        gene_dispersions = self.varm["genewise_dispersions"]
        use_for_mean = gene_dispersions > 10*self.min_disp
        mean_disp = scipy.stats.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)
```

Discussion:
scottgigante-immunai: @BorisMuzellec this is failing because IPython no longer supports Python3.8, unrelated to my PR. Leaving it to you to write tests to compare performance here to the R version.
BorisMuzellec: Thanks for this PR @scottgigante-immunai! I'll try to find some time this week to fix the IPython issue and write some R / Python tests.

Patch:
diff --git a/README.md b/README.md
index 8c70d6f..825f61b 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,9 @@ In case there is a feature you would particularly like to be implemented, feel f
 
 ## Installation
 
-`PyDESeq2` can be installed from PyPI:
+### PyPI
+
+`PyDESeq2` can be installed from PyPI using `pip`:
 
 `pip install pydeseq2`
 
@@ -46,6 +48,12 @@ conda install pip
 pip install pydeseq2
 ```
 
+### Bioconda
+
+`PyDESeq2` can also be installed from Bioconda with `conda`:
+
+`conda install -c bioconda pydeseq2`
+
 If you're interested in contributing or want access to the development version, please see the [contributing](#contributing) section.
 
 ### Requirements
@@ -97,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/docs/source/index.rst b/docs/source/index.rst
index e43391b..28d826c 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -15,7 +15,7 @@ As PyDESeq2 is a re-implementation of `DESeq2 <https://bioconductor.org/packages
 from scratch, you may experience some differences in terms of retrieved values or available features.
 
 Currently, available features broadly correspond to the default settings of DESeq2 (v1.34.0) for single-factor and
-paired multi-factor analysis (with bi-level factors), with an
+paired multi-factor analysis (with n-level categorical factors), with an
 optional `apeGLM <https://academic.oup.com/bioinformatics/article/35/12/2084/5159452>`_ LFC shrinkage step
 :cite:p:`zhu2019heavy`. We plan to implement more in the near future. In case there is a feature you would particularly
 like to be implemented, feel free to open an issue on GitHub.
diff --git a/pydeseq2/dds.py b/pydeseq2/dds.py
index 7ae1a35..887f135 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
@@ -226,7 +227,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 +242,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 +266,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.
diff --git a/pydeseq2/preprocessing.py b/pydeseq2/preprocessing.py
index 397f28c..2fa95ca 100644
--- a/pydeseq2/preprocessing.py
+++ b/pydeseq2/preprocessing.py
@@ -1,10 +1,13 @@
 from typing import Tuple
+from typing import Union
 
 import numpy as np
 import pandas as pd
 
 
-def deseq2_norm(counts: pd.DataFrame) -> Tuple[pd.DataFrame, pd.DataFrame]:
+def deseq2_norm(
+    counts: Union[pd.DataFrame, np.ndarray]
+) -> Tuple[Union[pd.DataFrame, np.ndarray], Union[pd.DataFrame, np.ndarray]]:
     """
     Return normalized counts and size_factors.
 
@@ -12,16 +15,16 @@ def deseq2_norm(counts: pd.DataFrame) -> Tuple[pd.DataFrame, pd.DataFrame]:
 
     Parameters
     ----------
-    counts : ndarray
+    counts : pandas.DataFrame or ndarray
             Raw counts. One column per gene, one row per sample.
 
     Returns
     -------
-    deseq2_counts : pandas.DataFrame
+    deseq2_counts : pandas.DataFrame or ndarray
         DESeq2 normalized counts.
         One column per gene, rows are indexed by sample barcodes.
 
-    size_factors : pandas.DataFrame
+    size_factors : pandas.DataFrame or ndarray
         DESeq2 normalization factors.
     """
 
@@ -32,7 +35,10 @@ def deseq2_norm(counts: pd.DataFrame) -> Tuple[pd.DataFrame, pd.DataFrame]:
     # Filter out genes with -∞ log means
     filtered_genes = ~np.isinf(logmeans)
     # Subtract filtered log means from log counts
-    log_ratios = log_counts[:, filtered_genes] - logmeans[filtered_genes]
+    if isinstance(log_counts, pd.DataFrame):
+        log_ratios = log_counts.loc[:, filtered_genes] - logmeans[filtered_genes]
+    else:
+        log_ratios = log_counts[:, filtered_genes] - logmeans[filtered_genes]
     # Compute sample-wise median of log ratios
     log_medians = np.median(log_ratios, axis=1)
     # Return raw counts divided by size factors (exponential of log ratios)


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_pydeseq2.py b/tests/test_pydeseq2.py
index 02cda60..0e6fcc2 100644
--- a/tests/test_pydeseq2.py
+++ b/tests/test_pydeseq2.py
@@ -8,6 +8,7 @@ import pandas as pd
 import tests
 from pydeseq2.dds import DeseqDataSet
 from pydeseq2.ds import DeseqStats
+from pydeseq2.preprocessing import deseq2_norm
 from pydeseq2.utils import load_example_data
 
 # Single-factor tests
@@ -448,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.
@@ -480,3 +514,34 @@ def test_ref_level():
             lambda x: 1 if x == "X" else 0 if x == "Y" else np.NaN
         )
     ).all()
+
+
+def test_deseq2_norm():
+    """Test that deseq2_norm() called on a pandas dataframe outputs the same results as
+    DeseqDataSet.fit_size_factors()
+    """
+    counts_df = load_example_data(
+        modality="raw_counts",
+        dataset="synthetic",
+        debug=False,
+    )
+
+    clinical_df = load_example_data(
+        modality="clinical",
+        dataset="synthetic",
+        debug=False,
+    )
+
+    # Fit size factors from DeseqDataSet
+    dds = DeseqDataSet(counts=counts_df, clinical=clinical_df)
+    dds.fit_size_factors()
+    s1 = dds.obsm["size_factors"]
+
+    # Fit size factors from counts directly
+    s2 = deseq2_norm(counts_df)[1]
+
+    np.testing.assert_almost_equal(
+        s1,
+        s2,
+        decimal=8,
+    )