refactored utils from de.testing into separate file

Author: davidsebfischer

diffxpy/api/test.py

@@ -1,2 +1,3 @@
-from diffxpy.testing.tests import design_matrix, coef_names, lrt, wald, t_test, rank_test, two_sample, pairwise, \
+from diffxpy.testing import lrt, wald, t_test, rank_test, two_sample, pairwise, \
     versus_rest, partition, continuous_1d
+from diffxpy.testing import design_matrix, coef_names

diffxpy/testing/__init__.py

@@ -0,0 +1,3 @@
+from .tests import lrt, wald, t_test, rank_test, two_sample, pairwise, \
+    versus_rest, partition, continuous_1d
+from .utils import design_matrix, coef_names

diffxpy/testing/det.py

@@ -7,7 +7,7 @@
 import xarray as xr
 import patsy
 
-from diffxpy.testing.tests import _split_X, t_test
+from .utils import split_X, dmat_unique
 
 try:
     import anndata
@@ -20,16 +20,7 @@
 from . import correction
 from diffxpy import pkg_constants
 
-logger = logging.getLogger(__name__)
-
-# Use this to suppress matrix subclass PendingDepreceationWarnings from numpy:
-np.warnings.filterwarnings("ignore")
-
-def _dmat_unique(dmat, sample_description):
-    dmat, idx = np.unique(dmat, axis=0, return_index=True)
-    sample_description = sample_description.iloc[idx].reset_index(drop=True)
-
-    return dmat, sample_description
+logger = logging.getLogger("diffxpy")
 
 
 class _Estimation(GeneralizedLinearModel, metaclass=abc.ABCMeta):
@@ -614,7 +605,7 @@ def _log_fold_change(self, factors: Union[Dict, Tuple, Set, List], base=np.e):
         dmat = self.full_estim.design_loc
 
         # make rows unique
-        dmat, sample_description = _dmat_unique(dmat, sample_description)
+        dmat, sample_description = dmat_unique(dmat, sample_description)
 
         # factors = factors.intersection(di.term_names)
 
@@ -628,7 +619,7 @@ def _log_fold_change(self, factors: Union[Dict, Tuple, Set, List], base=np.e):
         dmat[:, neg_sel] = 0
 
         # make the design matrix + sample description unique again
-        dmat, sample_description = _dmat_unique(dmat, sample_description)
+        dmat, sample_description = dmat_unique(dmat, sample_description)
 
         locations = self.full_estim.inverse_link_loc(dmat.dot(self.full_estim.par_link_loc))
         locations = np.log(locations) / np.log(base)
@@ -696,7 +687,7 @@ def locations(self):
         sample_description = self.sample_description[[f.name() for f in di.factor_infos]]
         dmat = self.full_estim.design_loc
 
-        dmat, sample_description = _dmat_unique(dmat, sample_description)
+        dmat, sample_description = dmat_unique(dmat, sample_description)
 
         retval = self.full_estim.inverse_link_loc(dmat.dot(self.full_estim.par_link_loc))
         retval = pd.DataFrame(retval, columns=self.full_estim.features)
@@ -718,7 +709,7 @@ def scales(self):
         sample_description = self.sample_description[[f.name() for f in di.factor_infos]]
         dmat = self.full_estim.design_scale
 
-        dmat, sample_description = _dmat_unique(dmat, sample_description)
+        dmat, sample_description = dmat_unique(dmat, sample_description)
 
         retval = self.full_estim.inverse_link_scale(dmat.doc(self.full_estim.par_link_scale))
         retval = pd.DataFrame(retval, columns=self.full_estim.features)
@@ -901,6 +892,7 @@ def summary(self, qval_thres=None, fc_upper_thres=None,
     def plot_vs_ttest(self, log10=False):
         import matplotlib.pyplot as plt
         import seaborn as sns
+        from .tests import t_test
 
         grouping = np.asarray(self.model_estim.design_loc[:, self.coef_loc_totest])
         ttest = t_test(
@@ -935,7 +927,7 @@ def __init__(self, data, grouping, gene_names, is_logged):
         self.grouping = grouping
         self._gene_names = np.asarray(gene_names)
 
-        x0, x1 = _split_X(data, grouping)
+        x0, x1 = split_X(data, grouping)
 
         # Only compute p-values for genes with non-zero observations and non-zero group-wise variance.
         mean_x0 = x0.mean(axis=0).astype(dtype=np.float)
@@ -1040,7 +1032,7 @@ def __init__(self, data, grouping, gene_names, is_logged):
         self.grouping = grouping
         self._gene_names = np.asarray(gene_names)
 
-        x0, x1 = _split_X(data, grouping)
+        x0, x1 = split_X(data, grouping)
 
         mean_x0 = x0.mean(axis=0).astype(dtype=np.float)
         mean_x1 = x1.mean(axis=0).astype(dtype=np.float)
@@ -1118,6 +1110,7 @@ def summary(self, qval_thres=None, fc_upper_thres=None,
     def plot_vs_ttest(self, log10=False):
         import matplotlib.pyplot as plt
         import seaborn as sns
+        from .tests import t_test
 
         grouping = self.grouping
         ttest = t_test(

diffxpy/testing/tests.py

@@ -1,5 +1,7 @@
 from typing import Union, List, Dict, Callable, Tuple
 
+import anndata
+import logging
 import numpy as np
 import pandas as pd
 import patsy
@@ -10,126 +12,17 @@
 from batchglm.xarray_sparse import SparseXArrayDataSet
 from diffxpy import pkg_constants
 from diffxpy.models.batch_bfgs.optim import Estim_BFGS
-from diffxpy.testing.det import anndata, logger, DifferentialExpressionTestLRT, DifferentialExpressionTestWald, \
+from .det import DifferentialExpressionTestLRT, DifferentialExpressionTestWald, \
     DifferentialExpressionTestTT, DifferentialExpressionTestRank, _DifferentialExpressionTestSingle, \
     DifferentialExpressionTestZTestLazy, DifferentialExpressionTestZTest, DifferentialExpressionTestPairwise, \
     DifferentialExpressionTestVsRest, _DifferentialExpressionTestMulti, DifferentialExpressionTestByPartition, \
     DifferentialExpressionTestWaldCont, DifferentialExpressionTestLRTCont
+from .utils import parse_gene_names, parse_data, parse_sample_description, parse_size_factors, parse_grouping
 
+logger = logging.getLogger("diffxpy")
 
-def _parse_gene_names(data, gene_names):
-    if gene_names is None:
-        if anndata is not None and (isinstance(data, anndata.AnnData) or isinstance(data, anndata.base.Raw)):
-            gene_names = data.var_names
-        elif isinstance(data, xr.DataArray):
-            gene_names = data["features"]
-        elif isinstance(data, xr.Dataset):
-            gene_names = data["features"]
-        else:
-            raise ValueError("Missing gene names")
-
-    return np.asarray(gene_names)
-
-
-def _parse_data(data, gene_names) -> xr.DataArray:
-    X = data_utils.xarray_from_data(data, dims=("observations", "features"))
-    if gene_names is not None:
-        X.coords["features"] = gene_names
-
-    return X
-
-
-def _parse_sample_description(data, sample_description=None) -> pd.DataFrame:
-    if sample_description is None:
-        if anndata is not None and isinstance(data, anndata.AnnData):
-            sample_description = data_utils.sample_description_from_anndata(
-                dataset=data,
-            )
-        elif isinstance(data, xr.Dataset):
-            sample_description = data_utils.sample_description_from_xarray(
-                dataset=data,
-                dim="observations",
-            )
-        else:
-            raise ValueError(
-                "Please specify `sample_description` or provide `data` as xarray.Dataset or anndata.AnnData " +
-                "with corresponding sample annotations"
-            )
-
-    if anndata is not None and isinstance(data, anndata.base.Raw):
-        # anndata.base.Raw does not have attribute shape.
-        assert data.X.shape[0] == sample_description.shape[0], \
-            "data matrix and sample description must contain same number of cells"
-    else:
-        assert data.shape[0] == sample_description.shape[0], \
-            "data matrix and sample description must contain same number of cells"
-    return sample_description
-
-
-def _parse_size_factors(size_factors, data):
-    if size_factors is not None:
-        if isinstance(size_factors, pd.core.series.Series):
-            size_factors = size_factors.values
-        assert size_factors.shape[0] == data.shape[0], "data matrix and size factors must contain same number of cells"
-    return size_factors
-
-
-def design_matrix(
-        data=None,
-        sample_description: pd.DataFrame = None,
-        formula: str = None,
-        dmat: pd.DataFrame = None
-) -> Union[patsy.design_info.DesignMatrix, xr.Dataset]:
-    """ Build design matrix for fit of generalized linear model.
-
-    This is necessary for wald tests and likelihood ratio tests.
-    This function only carries through formatting if dmat is directly supplied.
-
-    :param data: input data
-    :param formula: model formula.
-    :param sample_description: optional pandas.DataFrame containing sample annotations
-    :param dmat: model design matrix
-    """
-    if data is None and sample_description is None and dmat is None:
-        raise ValueError("Supply either data or sample_description or dmat.")
-    if dmat is None and formula is None:
-        raise ValueError("Supply either dmat or formula.")
-
-    if dmat is None:
-        sample_description = _parse_sample_description(data, sample_description)
-        dmat = data_utils.design_matrix(sample_description=sample_description, formula=formula)
-
-        return dmat
-    else:
-        ar = xr.DataArray(dmat, dims=("observations", "design_params"))
-        ar.coords["design_params"] = dmat.columns
-
-        ds = xr.Dataset({
-            "design": ar,
-        })
-
-        return ds
-
-
-def coef_names(
-        data=None,
-        sample_description: pd.DataFrame = None,
-        formula: str = None,
-        dmat: pd.DataFrame = None
-) -> list:
-    """ Output coefficient names of model only.
-
-    :param data: input data
-    :param formula: model formula.
-    :param sample_description: optional pandas.DataFrame containing sample annotations
-    :param dmat: model design matrix
-    """
-    return design_matrix(
-        data=data,
-        sample_description=sample_description,
-        formula=formula,
-        dmat=dmat
-    ).design_info.column_names
+# Use this to suppress matrix subclass PendingDepreceationWarnings from numpy:
+np.warnings.filterwarnings("ignore")
 
 
 def _fit(
@@ -408,10 +301,10 @@ def lrt(
     if isinstance(as_numeric, str):
         as_numeric = [as_numeric]
 
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
-    sample_description = _parse_sample_description(data, sample_description)
-    size_factors = _parse_size_factors(size_factors=size_factors, data=X)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
+    sample_description = parse_sample_description(data, sample_description)
+    size_factors = parse_size_factors(size_factors=size_factors, data=X)
 
     full_design_loc = data_utils.design_matrix(
         sample_description=sample_description,
@@ -623,11 +516,11 @@ def wald(
         as_numeric = [as_numeric]
 
     # # Parse input data formats:
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
     if dmat_loc is None and dmat_scale is None:
-        sample_description = _parse_sample_description(data, sample_description)
-    size_factors = _parse_size_factors(size_factors=size_factors, data=X)
+        sample_description = parse_sample_description(data, sample_description)
+    size_factors = parse_size_factors(size_factors=size_factors, data=X)
 
     if dmat_loc is None:
         design_loc = data_utils.design_matrix(
@@ -714,20 +607,6 @@ def wald(
     return de_test
 
 
-def _parse_grouping(data, sample_description, grouping):
-    if isinstance(grouping, str):
-        sample_description = _parse_sample_description(data, sample_description)
-        grouping = sample_description[grouping]
-    return np.squeeze(np.asarray(grouping))
-
-
-def _split_X(data, grouping):
-    groups = np.unique(grouping)
-    x0 = data[np.where(grouping == groups[0])[0]]
-    x1 = data[np.where(grouping == groups[1])[0]]
-    return x0, x1
-
-
 def t_test(
         data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix],
         grouping,
@@ -752,11 +631,11 @@ def t_test(
         Whether data is already logged. If True, log-fold changes are computed as fold changes on this data.
         If False, log-fold changes are computed as log-fold changes on this data.
     """
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
     if isinstance(X, SparseXArrayDataSet):
         X = X.X
-    grouping = _parse_grouping(data, sample_description, grouping)
+    grouping = parse_grouping(data, sample_description, grouping)
 
     de_test = DifferentialExpressionTestTT(
         data=X.astype(dtype),
@@ -792,11 +671,11 @@ def rank_test(
         Whether data is already logged. If True, log-fold changes are computed as fold changes on this data.
         If False, log-fold changes are computed as log-fold changes on this data.
     """
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
     if isinstance(X, SparseXArrayDataSet):
         X = X.X
-    grouping = _parse_grouping(data, sample_description, grouping)
+    grouping = parse_grouping(data, sample_description, grouping)
 
     de_test = DifferentialExpressionTestRank(
         data=X.astype(dtype),
@@ -910,9 +789,9 @@ def two_sample(
         raise ValueError('base.two_sample(): Do not specify `noise_model` if using test t-test or wilcoxon: ' +
                          'The t-test is based on a gaussian noise model and wilcoxon is model free.')
 
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
-    grouping = _parse_grouping(data, sample_description, grouping)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
+    grouping = parse_grouping(data, sample_description, grouping)
     sample_description = pd.DataFrame({"grouping": grouping})
 
     groups = np.unique(grouping)
@@ -1116,10 +995,10 @@ def pairwise(
 
     # Do not store all models but only p-value and q-value matrix:
     # genes x groups x groups
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
-    sample_description = _parse_sample_description(data, sample_description)
-    grouping = _parse_grouping(data, sample_description, grouping)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
+    sample_description = parse_sample_description(data, sample_description)
+    grouping = parse_grouping(data, sample_description, grouping)
     sample_description = pd.DataFrame({"grouping": grouping})
 
     if test.lower() == 'z-test' or test.lower() == 'z_test' or test.lower() == 'ztest':
@@ -1324,10 +1203,10 @@ def versus_rest(
 
     # Do not store all models but only p-value and q-value matrix:
     # genes x groups
-    gene_names = _parse_gene_names(data, gene_names)
-    X = _parse_data(data, gene_names)
-    sample_description = _parse_sample_description(data, sample_description)
-    grouping = _parse_grouping(data, sample_description, grouping)
+    gene_names = parse_gene_names(data, gene_names)
+    X = parse_data(data, gene_names)
+    sample_description = parse_sample_description(data, sample_description)
+    grouping = parse_grouping(data, sample_description, grouping)
     sample_description = pd.DataFrame({"grouping": grouping})
 
     groups = np.unique(grouping)
@@ -1427,10 +1306,10 @@ def __init__(
         :param gene_names: optional list/array of gene names which will be used if `data` does not implicitly store these
         :param sample_description: optional pandas.DataFrame containing sample annotations
         """
-        self.X = _parse_data(data, gene_names)
-        self.gene_names = _parse_gene_names(data, gene_names)
-        self.sample_description = _parse_sample_description(data, sample_description)
-        self.partition = _parse_grouping(data, sample_description, partition)
+        self.X = parse_data(data, gene_names)
+        self.gene_names = parse_gene_names(data, gene_names)
+        self.sample_description = parse_sample_description(data, sample_description)
+        self.partition = parse_grouping(data, sample_description, partition)
         self.partitions = np.unique(self.partition)
         self.partition_idx = [np.where(self.partition == x)[0] for x in self.partitions]
 
@@ -1874,9 +1753,9 @@ def continuous_1d(
     if isinstance(as_numeric, tuple):
         as_numeric = list(as_numeric)
 
-    X = _parse_data(data, gene_names)
-    gene_names = _parse_gene_names(data, gene_names)
-    sample_description = _parse_sample_description(data, sample_description)
+    X = parse_data(data, gene_names)
+    gene_names = parse_gene_names(data, gene_names)
+    sample_description = parse_sample_description(data, sample_description)
 
     # Check that continuous factor is contained in sample description
     if continuous not in sample_description.columns: