Add permutation test

docs/usage/usage.md

@@ -136,6 +136,7 @@ Pertpy provides utilities to conduct differential gene expression tests through
     tools.EdgeR
     tools.WilcoxonTest
     tools.TTest
+    tools.PermutationTest
     tools.Statsmodels
 ```
 
@@ -563,33 +564,33 @@ including cell line annotation, bulk RNA and protein expression data.
 
 Available databases for cell line metadata:
 
-- [The Cancer Dependency Map Project at Broad](https://depmap.org/portal/)
-- [The Cancer Dependency Map Project at Sanger](https://depmap.sanger.ac.uk/)
-- [Genomics of Drug Sensitivity in Cancer (GDSC)](https://www.cancerrxgene.org/)
+-   [The Cancer Dependency Map Project at Broad](https://depmap.org/portal/)
+-   [The Cancer Dependency Map Project at Sanger](https://depmap.sanger.ac.uk/)
+-   [Genomics of Drug Sensitivity in Cancer (GDSC)](https://www.cancerrxgene.org/)
 
 ### Compound
 
 The Compound module enables the retrieval of various types of information related to compounds of interest, including the most common synonym, pubchemID and canonical SMILES.
 
 Available databases for compound metadata:
 
-- [PubChem](https://pubchem.ncbi.nlm.nih.gov/)
+-   [PubChem](https://pubchem.ncbi.nlm.nih.gov/)
 
 ### Mechanism of Action
 
 This module aims to retrieve metadata of mechanism of action studies related to perturbagens of interest, depending on the molecular targets.
 
 Available databases for mechanism of action metadata:
 
-- [CLUE](https://clue.io/)
+-   [CLUE](https://clue.io/)
 
 ### Drug
 
 This module allows for the retrieval of Drug target information.
 
 Available databases for drug metadata:
 
-- [chembl](https://www.ebi.ac.uk/chembl/)
+-   [chembl](https://www.ebi.ac.uk/chembl/)
 
 ```{eval-rst}
 .. autosummary::

pertpy/tools/__init__.py

@@ -47,6 +47,7 @@ def __init__(self, *args, **kwargs):
 
 DE_EXTRAS = ["formulaic", "pydeseq2"]
 EdgeR = lazy_import("pertpy.tools._differential_gene_expression", "EdgeR", DE_EXTRAS)  # edgeR will be imported via rpy2
+PermutationTest = lazy_import("pertpy.tools._differential_gene_expression", "PermutationTest", DE_EXTRAS)
 PyDESeq2 = lazy_import("pertpy.tools._differential_gene_expression", "PyDESeq2", DE_EXTRAS)
 Statsmodels = lazy_import("pertpy.tools._differential_gene_expression", "Statsmodels", DE_EXTRAS + ["statsmodels"])
 TTest = lazy_import("pertpy.tools._differential_gene_expression", "TTest", DE_EXTRAS)
@@ -62,6 +63,7 @@ def __init__(self, *args, **kwargs):
     "PyDESeq2",
     "WilcoxonTest",
     "TTest",
+    "PermutationTest",
     "Statsmodels",
     "DistanceTest",
     "Distance",

pertpy/tools/_differential_gene_expression/__init__.py

@@ -2,7 +2,7 @@
 from ._dge_comparison import DGEEVAL
 from ._edger import EdgeR
 from ._pydeseq2 import PyDESeq2
-from ._simple_tests import SimpleComparisonBase, TTest, WilcoxonTest
+from ._simple_tests import PermutationTest, SimpleComparisonBase, TTest, WilcoxonTest
 from ._statsmodels import Statsmodels
 
 __all__ = [
@@ -14,6 +14,7 @@
     "SimpleComparisonBase",
     "WilcoxonTest",
     "TTest",
+    "PermutationTest",
 ]
 
-AVAILABLE_METHODS = [Statsmodels, EdgeR, PyDESeq2, WilcoxonTest, TTest]
+AVAILABLE_METHODS = [Statsmodels, EdgeR, PyDESeq2, WilcoxonTest, TTest, PermutationTest]

pertpy/tools/_differential_gene_expression/_pydeseq2.py

@@ -42,7 +42,7 @@ def fit(self, **kwargs) -> pd.DataFrame:
             **kwargs: Keyword arguments specific to DeseqDataSet(), except for `n_cpus` which will use all available CPUs minus one if the argument is not passed.
         """
         try:
-            usable_cpus = len(os.sched_getaffinity(0))
+            usable_cpus = len(os.sched_getaffinity(0))  # type: ignore # os.sched_getaffinity is not available on Windows and macOS
         except AttributeError:
             usable_cpus = os.cpu_count()
 

pertpy/tools/_differential_gene_expression/_simple_tests.py

@@ -2,14 +2,16 @@
 
 import warnings
 from abc import abstractmethod
-from collections.abc import Mapping, Sequence
+from collections.abc import Callable, Mapping, Sequence
 from types import MappingProxyType
 
 import numpy as np
 import pandas as pd
 import scipy.stats
 import statsmodels
 from anndata import AnnData
+from joblib import Parallel, delayed
+from lamin_utils import logger
 from pandas.core.api import DataFrame as DataFrame
 from scipy.sparse import diags, issparse
 from tqdm.auto import tqdm
@@ -94,9 +96,28 @@ def compare_groups(
         paired_by: str | None = None,
         mask: str | None = None,
         layer: str | None = None,
+        n_permutations: int = 1000,
+        permutation_test: type["SimpleComparisonBase"] | None = None,
         fit_kwargs: Mapping = MappingProxyType({}),
         test_kwargs: Mapping = MappingProxyType({}),
+        n_jobs: int = -1,
     ) -> DataFrame:
+        """Perform a comparison between groups.
+
+        Args:
+            adata (AnnData): Data with observations to compare.
+            column (str): Column in `adata.obs` that contains the groups to compare.
+            baseline (str): Reference group.
+            groups_to_compare (str | Sequence[str]): Groups to compare against the baseline. If None, all other groups are compared.
+            paired_by (str | None): Column in `adata.obs` to use for pairing. If None, an unpaired test is performed.
+            mask (str | None): Mask to apply to the data.
+            layer (str | None): Layer to use for the comparison.
+            n_permutations (int): Number of permutations to perform if a permutation test is used.
+            permutation_test (type[SimpleComparisonBase] | None): Test to use after permutation if a permutation test is used.
+            fit_kwargs (Mapping): Not used for simple tests.
+            test_kwargs (Mapping): Additional kwargs passed to the test function.
+            n_jobs (int): Number of parallel jobs to use.
+        """
         if len(fit_kwargs):
             warnings.warn("fit_kwargs not used for simple tests.", UserWarning, stacklevel=2)
         paired = paired_by is not None
@@ -127,13 +148,24 @@ def _get_idx(column, value):
 
         res_dfs = []
         baseline_idx = _get_idx(column, baseline)
-        for group_to_compare in groups_to_compare:
-            comparison_idx = _get_idx(column, group_to_compare)
-            res_dfs.append(
-                model._compare_single_group(baseline_idx, comparison_idx, paired=paired, **test_kwargs).assign(
-                    comparison=f"{group_to_compare}_vs_{baseline if baseline is not None else 'rest'}"
-                )
+
+        if permutation_test:
+            test_kwargs = dict(test_kwargs)
+            test_kwargs.update({"test": permutation_test, "n_permutations": n_permutations})
+        elif permutation_test is None and cls.__name__ == "PermutationTest":
+            logger.warning("No permutation test specified. Using WilcoxonTest as default.")
+
+        comparison_indices = [_get_idx(column, group_to_compare) for group_to_compare in groups_to_compare]
+        res_dfs = Parallel(n_jobs=n_jobs)(
+            delayed(model._compare_single_group)(baseline_idx, comparison_idx, paired=paired, **test_kwargs)
+            for comparison_idx in comparison_indices
+        )
+        res_dfs = [
+            df.assign(
+                comparison=f"{group_to_compare}_vs_{baseline if baseline is not None else 'rest'}",
             )
+            for df, group_to_compare in zip(res_dfs, groups_to_compare, strict=False)
+        ]
         return fdr_correction(pd.concat(res_dfs))
 
 
@@ -144,19 +176,100 @@ class WilcoxonTest(SimpleComparisonBase):
     """
 
     @staticmethod
-    def _test(x0: np.ndarray, x1: np.ndarray, paired: bool, **kwargs) -> float:
+    def _test(x0: np.ndarray, x1: np.ndarray, paired: bool, return_attribute: str = "pvalue", **kwargs) -> float:
         if paired:
-            return scipy.stats.wilcoxon(x0, x1, **kwargs).pvalue
+            return scipy.stats.wilcoxon(x0, x1, **kwargs).__getattribute__(return_attribute)
         else:
-            return scipy.stats.mannwhitneyu(x0, x1, **kwargs).pvalue
+            return scipy.stats.mannwhitneyu(x0, x1, **kwargs).__getattribute__(return_attribute)
 
 
 class TTest(SimpleComparisonBase):
-    """Perform a unpaired or paired T-test"""
+    """Perform a unpaired or paired T-test."""
 
     @staticmethod
-    def _test(x0: np.ndarray, x1: np.ndarray, paired: bool, **kwargs) -> float:
+    def _test(x0: np.ndarray, x1: np.ndarray, paired: bool, return_attribute: str = "pvalue", **kwargs) -> float:
         if paired:
-            return scipy.stats.ttest_rel(x0, x1, **kwargs).pvalue
+            return scipy.stats.ttest_rel(x0, x1, **kwargs).__getattribute__(return_attribute)
         else:
-            return scipy.stats.ttest_ind(x0, x1, **kwargs).pvalue
+            return scipy.stats.ttest_ind(x0, x1, **kwargs).__getattribute__(return_attribute)
+
+
+class PermutationTest(SimpleComparisonBase):
+    """Perform a permutation test.
+
+    The permutation test relies on another test (e.g. WilcoxonTest) to perform the actual comparison
+    based on permuted data. The p-value is then calculated based on the distribution of the test
+    statistic under the null hypothesis.
+
+    For paired tests, each paired observation is permuted together and distributed randomly between
+    the two groups. For unpaired tests, all observations are permuted independently.
+
+    The null hypothesis for the unpaired test is that all observations come from the same underlying
+    distribution and have been randomly assigned to one of the samples.
+
+    The null hypothesis for the paired permutation test is that the observations within each pair are
+    drawn from the same underlying distribution and that their assignment to a sample is random.
+    """
+
+    @staticmethod
+    def _test(
+        x0: np.ndarray,
+        x1: np.ndarray,
+        paired: bool,
+        test: type["SimpleComparisonBase"] | Callable = WilcoxonTest,
+        n_permutations: int = 1000,
+        return_attribute: str = "pvalue",
+        **kwargs,
+    ) -> float:
+        """Perform a permutation test.
+
+        This function relies on another test (e.g. WilcoxonTest) to generate a test statistic for each permutation.
+
+        .. code-block:: python
+            from pertpy.tools import PermutationTest, WilcoxonTest
+
+            # Using rank-sum statistic
+            p_value = PermutationTest._test(x0, x1, paired=True, test=WilcoxonTest, n_permutations=1000, rng=0)
+
+
+            # Using a custom test statistic
+            def compare_means(x0, x1, paired):
+                # paired logic not implemented here
+                return np.mean(x1) - np.mean(x0)
+
+
+            p_value = PermutationTest._test(x0, x1, paired=False, test=compare_means, n_permutations=1000, rng=0)
+
+        Args:
+            x0: Array with baseline values.
+            x1: Array with values to compare.
+            paired: Whether to perform a paired test
+            test: The class or function to generate the test statistic from permuted data.
+            n_permutations: Number of permutations to perform.
+            return_attribute: Attribute to return from the test statistic.
+            **kwargs: kwargs passed to the permutation test function, not the test function after permutation.
+        """
+        if test is PermutationTest:
+            raise ValueError(
+                "The `test` argument cannot be `PermutationTest`. Use a base test like `WilcoxonTest` or `TTest`."
+            )
+
+        def call_test(data_baseline, data_comparison, axis: int | None = None, **kwargs):
+            """Perform the actual test."""
+            # Setting the axis allows the operation to be vectorized
+            if axis is not None:
+                kwargs.update({"axis": axis})
+
+            if not hasattr(test, "_test"):
+                return test(data_baseline, data_comparison, paired, **kwargs)
+
+            return test._test(data_baseline, data_comparison, paired, return_attribute="statistic", **kwargs)
+
+        return scipy.stats.permutation_test(
+            [x0, x1],
+            statistic=call_test,
+            n_resamples=n_permutations,
+            permutation_type=("samples" if paired else "independent"),
+            vectorized=hasattr(test, "_test"),
+            **kwargs,
+        ).__getattribute__(return_attribute)

tests/tools/_differential_gene_expression/test_simple_tests.py

@@ -2,7 +2,7 @@
 import pandas as pd
 import pytest
 from pandas.core.api import DataFrame as DataFrame
-from pertpy.tools._differential_gene_expression import SimpleComparisonBase, TTest, WilcoxonTest
+from pertpy.tools._differential_gene_expression import PermutationTest, SimpleComparisonBase, TTest, WilcoxonTest
 
 
 @pytest.mark.parametrize(
@@ -61,6 +61,45 @@ def test_t(test_adata_minimal, paired_by, expected):
         assert actual[gene] == pytest.approx(expected[gene], abs=0.02)
 
 
+@pytest.mark.parametrize(
+    "paired_by,expected",
+    [
+        pytest.param(
+            None,
+            {"gene1": {"p_value": 2.13e-26, "log_fc": -5.14}, "gene2": {"p_value": 0.96, "log_fc": -0.016}},
+            id="unpaired",
+        ),
+        pytest.param(
+            "pairing",
+            {"gene1": {"p_value": 1.63e-26, "log_fc": -5.14}, "gene2": {"p_value": 0.85, "log_fc": -0.016}},
+            id="paired",
+        ),
+    ],
+)
+def test_permutation(test_adata_minimal, paired_by, expected):
+    """Test that permutation test gives the correct values.
+
+    Reference values have been computed in R using wilcox.test
+    """
+    for permutation_test in [TTest, WilcoxonTest]:
+        res_df = PermutationTest.compare_groups(
+            adata=test_adata_minimal,
+            column="condition",
+            baseline="A",
+            groups_to_compare="B",
+            paired_by=paired_by,
+            n_permutations=200,
+            permutation_test=permutation_test,
+            test_kwargs={"rng": 0},
+        )
+        assert isinstance(res_df, DataFrame), "PermutationTest.compare_groups should return a DataFrame"
+        actual = res_df.loc[:, ["variable", "p_value", "log_fc"]].set_index("variable").to_dict(orient="index")
+        for gene in expected:
+            assert (actual[gene]["p_value"] < 0.05) == (expected[gene]["p_value"] < 0.05)
+            if actual[gene]["p_value"] < 0.05:
+                assert actual[gene] == pytest.approx(expected[gene], abs=0.02)
+
+
 @pytest.mark.parametrize("seed", range(10))
 def test_simple_comparison_pairing(test_adata_minimal, seed):
     """Test that paired samples are properly matched in a paired test"""