refactored single test unit tests to ease phasing in of new noise models

davidsebfischer · davidsebfischer · commit b851187d9a08 · 2019-04-03T10:12:54.000+02:00
diff --git a/diffxpy/unit_test/test_single_de.py b/diffxpy/unit_test/test_single_de.py
@@ -0,0 +1,268 @@
+import unittest
+import logging
+import numpy as np
+import pandas as pd
+import scipy.stats as stats
+
+from batchglm.api.models.glm_nb import Simulator
+import diffxpy.api as de
+
+
+class _TestSingleDE:
+
+    def _prepare_data(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        num_non_de = n_genes // 2
+        sim = Simulator(num_observations=n_cells, num_features=n_genes)
+        sim.generate_sample_description(num_batches=0, num_conditions=2)
+        sim.generate_params(
+            rand_fn_ave=lambda shape: np.random.poisson(500, shape) + 1,
+            rand_fn=lambda shape: np.abs(np.random.uniform(1, 0.5, shape))
+        )
+        sim.params["a_var"][1, :num_non_de] = 0
+        sim.params["b_var"][1, :num_non_de] = 0
+        sim.params["isDE"] = ("features",), np.arange(n_genes) >= num_non_de
+        sim.generate_data()
+
+        return sim
+
+    def _eval(self, sim, test):
+        idx_de = np.where(sim.params["isDE"] == True)[0]
+        idx_nonde = np.where(sim.params["isDE"] == False)[0]
+
+        frac_de_of_non_de = np.sum(test.qval[idx_nonde] < 0.05) / len(idx_nonde)
+        frac_de_of_de = np.sum(test.qval[idx_de] < 0.05) / len(idx_de)
+
+        logging.getLogger("diffxpy").info(
+            'fraction of non-DE genes with q-value < 0.05: %.1f%%' %
+            float(100 * frac_de_of_non_de)
+        )
+        logging.getLogger("diffxpy").info(
+            'fraction of DE genes with q-value < 0.05: %.1f%%' %
+            float(100 * frac_de_of_de)
+        )
+        assert frac_de_of_non_de <= 0.1, "too many false-positives"
+        assert frac_de_of_de >= 0.5, "too many false-negatives"
+
+        return sim
+
+    def _test_rank_de(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
+
+        sim = self._prepare_data(n_cells=n_cells, n_genes=n_genes)
+
+        test = de.test.rank_test(
+            data=sim.X,
+            grouping="condition",
+            sample_description=sim.sample_description,
+            dtype="float64"
+        )
+
+        self._eval(sim=sim, test=test)
+
+        return True
+
+    def _test_t_test_de(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
+
+        sim = self._prepare_data(n_cells=n_cells, n_genes=n_genes)
+
+        test = de.test.t_test(
+            data=sim.X,
+            grouping="condition",
+            sample_description=sim.sample_description,
+            dtype="float64"
+        )
+
+        self._eval(sim=sim, test=test)
+
+        return True
+
+    def _test_wald_de(
+            self,
+            n_cells: int,
+            n_genes: int,
+            noise_model: str
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        :param noise_model: Noise model to use for data fitting.
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
+
+        sim = self._prepare_data(n_cells=n_cells, n_genes=n_genes)
+
+        test = de.test.wald(
+            data=sim.X,
+            factor_loc_totest="condition",
+            formula_loc="~ 1 + condition",
+            sample_description=sim.sample_description,
+            noise_model=noise_model,
+            training_strategy="DEFAULT",
+            dtype="float64"
+        )
+
+        self._eval(sim=sim, test=test)
+
+        return True
+
+    def _test_lrt_de(
+            self,
+            n_cells: int,
+            n_genes: int,
+            noise_model: str
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        :param noise_model: Noise model to use for data fitting.
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
+
+        sim = self._prepare_data(n_cells=n_cells, n_genes=n_genes)
+
+        test = de.test.lrt(
+            data=sim.X,
+            full_formula_loc="~ 1 + condition",
+            full_formula_scale="~ 1",
+            reduced_formula_loc="~ 1",
+            reduced_formula_scale="~ 1",
+            sample_description=sim.sample_description,
+            noise_model=noise_model,
+            training_strategy="DEFAULT",
+            dtype="float64"
+        )
+
+        self._eval(sim=sim, test=test)
+
+        return True
+
+
+class TestSingleDE_STANDARD(_TestSingleDE, unittest.TestCase):
+    """
+    Noise model-independent tests unit tests that tests false positive and false negative rates.
+    """
+
+    def test_ttest_de(
+            self,
+            n_cells: int = 2000,
+            n_genes: int = 200
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        return self._test_t_test_de(
+            n_cells=n_cells,
+            n_genes=n_genes
+        )
+
+    def test_rank_de(
+            self,
+            n_cells: int = 2000,
+            n_genes: int = 200
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        return self._test_rank_de(
+            n_cells=n_cells,
+            n_genes=n_genes
+        )
+
+
+class TestSingleDE_NB(_TestSingleDE, unittest.TestCase):
+    """
+    Negative binomial noise model unit tests that tests false positive and false negative rates.
+    """
+
+    def test_wald_de_nb(
+            self,
+            n_cells: int = 2000,
+            n_genes: int = 200
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        return self._test_wald_de(
+            n_cells=n_cells,
+            n_genes=n_genes,
+            noise_model="nb"
+        )
+
+    def test_lrt_de_nb(
+            self,
+            n_cells: int = 2000,
+            n_genes: int = 200
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        return self._test_lrt_de(
+            n_cells=n_cells,
+            n_genes=n_genes,
+            noise_model="nb"
+        )
+
+
+class TestSingleDE_NORM(_TestSingleDE, unittest.TestCase):
+    """
+    Normal noise model unit tests that tests false positive and false negative rates.
+    """
+
+    def test_wald_de_norm(
+            self,
+            n_cells: int = 2000,
+            n_genes: int = 200
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        return self._test_wald_de(
+            n_cells=n_cells,
+            n_genes=n_genes,
+            noise_model="norm"
+        )
+
+    def test_lrt_de_norm(
+            self,
+            n_cells: int = 2000,
+            n_genes: int = 200
+    ):
+        """
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        return self._test_lrt_de(
+            n_cells=n_cells,
+            n_genes=n_genes,
+            noise_model="norm"
+        )
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/diffxpy/unit_test/test_single_external_libs.py b/diffxpy/unit_test/test_single_external_libs.py
@@ -0,0 +1,117 @@
+import unittest
+import logging
+import numpy as np
+import pandas as pd
+import scipy.stats as stats
+
+from batchglm.api.models.glm_nb import Simulator
+import diffxpy.api as de
+
+
+class TestSingleExternalLibs(unittest.TestCase):
+
+    def _prepare_data(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        sim = Simulator(num_observations=n_cells, num_features=n_genes)
+        sim.generate_sample_description(num_batches=0, num_conditions=2)
+        sim.generate_params(
+            rand_fn_ave=lambda shape: np.random.poisson(500, shape) + 1,
+            rand_fn=lambda shape: np.abs(np.random.uniform(1, 0.5, shape))
+        )
+        sim.generate_data()
+
+        return sim
+
+    def _eval(self, test, ref_pvals):
+        test_pval = test.pval
+        pval_dev = np.abs(test_pval - ref_pvals)
+        log_pval_dev = np.abs(np.log(test_pval+1e-200) - np.log(ref_pvals+1e-200))
+        max_dev = np.max(pval_dev)
+        max_log_dev = np.max(log_pval_dev)
+        mean_dev = np.mean(log_pval_dev)
+        logging.getLogger("diffxpy").info(
+            'maximum absolute p-value deviation: %f' %
+            float(max_dev)
+        )
+        logging.getLogger("diffxpy").info(
+            'maximum absolute log p-value deviation: %f' %
+            float(max_log_dev)
+        )
+        logging.getLogger("diffxpy").info(
+            'mean absolute log p-value deviation: %f' %
+            float(mean_dev)
+        )
+        assert max_dev < 1e-3, "maximum deviation too large"
+        assert max_log_dev < 1e-1, "maximum deviation in log space too large"
+
+    def test_t_test_ref(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+        Test if de.test.t_test() generates the same p-value distribution as scipy t-test.
+
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.INFO)
+
+        sim = self._prepare_data(n_cells=n_cells, n_genes=n_genes)
+
+        test = de.test.t_test(
+            data=sim.X,
+            grouping="condition",
+            sample_description=sim.sample_description,
+            dtype="float64"
+        )
+
+        # Run scipy t-tests as a reference.
+        conds = np.unique(sim.sample_description["condition"].values)
+        ind_a = np.where(sim.sample_description["condition"] == conds[0])[0]
+        ind_b = np.where(sim.sample_description["condition"] == conds[1])[0]
+        scipy_pvals = stats.ttest_ind(a=sim.X[ind_a, :], b=sim.X[ind_b, :], axis=0, equal_var=False).pvalue
+
+        self._eval(test=test, ref_pvals=scipy_pvals)
+
+        return True
+
+    def test_rank_ref(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+        Test if de.test.rank_test() generates the same p-value distribution as scipy t-test.
+
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.INFO)
+
+        sim = self._prepare_data(n_cells=n_cells, n_genes=n_genes)
+
+        test = de.test.rank_test(
+            data=sim.X,
+            grouping="condition",
+            sample_description=sim.sample_description,
+            dtype="float64"
+        )
+
+        # Run scipy t-tests as a reference.
+        conds = np.unique(sim.sample_description["condition"].values)
+        ind_a = np.where(sim.sample_description["condition"] == conds[0])[0]
+        ind_b = np.where(sim.sample_description["condition"] == conds[1])[0]
+        scipy_pvals = np.array([
+            stats.mannwhitneyu(x=sim.X[ind_a, i], y=sim.X[ind_b, i],
+                               use_continuity=True, alternative="two-sided").pvalue
+            for i in range(sim.X.shape[1])
+            ])
+
+        self._eval(test=test, ref_pvals=scipy_pvals)
+
+        return True
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/diffxpy/unit_test/test_single_null.py b/diffxpy/unit_test/test_single_null.py
