diff --git a/LICENSE b/LICENSE index 0d8d25e..db8e2b4 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ BSD 3-Clause License -Copyright (c) 2018, David S. Fischer, Florian R. Hölzlwimmer. +Copyright (c) 2018, David S. Fischer, Florian R. Hölzlwimmer, Theis Lab (theislab). All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/NOTICE b/NOTICE new file mode 100644 index 0000000..e3e6df7 --- /dev/null +++ b/NOTICE @@ -0,0 +1,2 @@ +The file ./docs/conf.py was adapted from ttps://github.com/theislab/scanpy/scanpy/conf.py and is licensed by the +scanpy project (F. Alexander Wolf, P. Angerer, Theis Lab) as described in the file. \ No newline at end of file diff --git a/diffxpy/__init__.py b/diffxpy/__init__.py index 2effa51..a2b32b8 100644 --- a/diffxpy/__init__.py +++ b/diffxpy/__init__.py @@ -4,3 +4,11 @@ del get_versions from .log_cfg import logger, unconfigure_logging, enable_logging + +__author__ = ', '.join([ + 'David Sebastian Fischer', + 'Florian Hölzlwimmer' +]) +__email__ = ', '.join([ + 'david.fischer@helmholtz-muenchen.de' +]) diff --git a/diffxpy/api/test.py b/diffxpy/api/test.py index e44a2c8..95ff885 100644 --- a/diffxpy/api/test.py +++ b/diffxpy/api/test.py @@ -1,3 +1,2 @@ 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 diff --git a/diffxpy/api/utils.py b/diffxpy/api/utils.py index c4acbec..fb653a7 100644 --- a/diffxpy/api/utils.py +++ b/diffxpy/api/utils.py @@ -1 +1,4 @@ -import batchglm.data as data_utils \ No newline at end of file +from diffxpy.testing.utils import constraint_matrix_from_string, constraint_matrix_from_dict, \ + constraint_system_from_star +from diffxpy.testing.utils import design_matrix, design_matrix_from_xarray, design_matrix_from_anndata +from diffxpy.testing.utils import view_coef_names, preview_coef_names diff --git a/diffxpy/enrichment/enrich.py b/diffxpy/enrichment/enrich.py index 812cb46..f1d4a1d 100644 --- a/diffxpy/enrichment/enrich.py +++ b/diffxpy/enrichment/enrich.py @@ -7,7 +7,6 @@ from ..testing import correction from ..testing.det import _DifferentialExpressionTest -logger = logging.getLogger(__name__) class RefSets: """ @@ -42,15 +41,19 @@ def clean(self, ids): def __init__(self, sets=None, fn=None, type='gmt'): if sets is not None: - self.load_sets(sets, type=type) - self._genes = np.sort(np.unique(np.concatenate([np.asarray(list(x.genes)) for x in self.sets]))) + if len(sets) > 0: + self.load_sets(sets, type=type) + self._genes = np.sort(np.unique(np.concatenate([np.asarray(list(x.genes)) for x in self.sets]))) + else: + self.sets = [] + self._genes = np.array([]) elif fn is not None: self.read_from_file(fn=fn, type=type) self._genes = np.sort(np.unique(np.concatenate([np.asarray(list(x.genes)) for x in self.sets]))) else: self.sets = [] self._genes = np.array([]) - self._ids = [x.id for x in self.sets] + self._ids = np.array([x.id for x in self.sets]) self._set_lens = np.array([x.len for x in self.sets]) self.genes_discarded = None @@ -113,7 +116,7 @@ def add(self, id: str, source: str, gene_ids: list): self.sets.append(self._Set(id=id, source=source, gene_ids=gene_ids)) # Update summary variables: self._genes = np.sort(np.unique(np.concatenate([np.asarray(list(x.genes)) for x in self.sets]))) - self._ids = [x.id for x in self.sets] + self._ids = np.array([x.id for x in self.sets]) self._set_lens = np.array([x.len for x in self.sets]) ## Processing functions. @@ -165,7 +168,7 @@ def get_set(self, id): """ Return the set with a given set identifier. """ - return self.sets[self._ids.index(id)] + return self.sets[self._ids.tolist().index(id)] ## Overlap functions. @@ -191,9 +194,9 @@ def overlap(self, enq_set: set, set_id=None): def test( ref: RefSets, det: Union[_DifferentialExpressionTest, None] = None, - pval: Union[np.array, None] = None, + scores: Union[np.array, None] = None, gene_ids: Union[list, None] = None, - de_threshold=0.05, + threshold=0.05, incl_all_zero=False, all_ids=None, clean_ref=False, @@ -205,38 +208,31 @@ def test( nice doc string and that the call to this is de.enrich.test which makes more sense to me than de.enrich.Enrich. - :param RefSets: - The annotated gene sets against which enrichment is tested. - :param DETest: - The differential expression results object which is tested + :param ref: The annotated gene sets against which enrichment is tested. + :param det: The differential expression results object which is tested for enrichment in the gene sets. - :param pval: - Alternative to DETest, vector of p-values for differential expression. - :param gene_ids: - If pval was supplied instead of DETest, use gene_ids to supply the + :param scores: Alternative to DETest, vector of scores (scalar per gene) which are then + used to discretize gene list. This can for example be corrected p-values from a differential expression + test, in that case the parameter threshold would be a significance threshold. + :param gene_ids: If pval was supplied instead of DETest, use gene_ids to supply the vector of gene identifiers (strings) that correspond to the p-values - which can be matched against the identifieres in the sets in RefSets. - :param de_threshold: - Significance threshold at which a differential test (a multiple-testing - corrected p-value) is called siginficant. T - :param incl_all_zero: - Wehther to include genes in gene universe which were all zero. - :param all_ids: - Set of all gene identifiers, this is used as the background set in the + which can be matched against the identifiers in the sets in RefSets. + :param threshold: Threshold of parameter scores at which a gene is included as a hit: In the case + of differential test p-values in scores, threshold is the significance threshold. + :param incl_all_zero: Wehther to include genes in gene universe which were all zero. + :param all_ids: Set of all gene identifiers, this is used as the background set in the hypergeometric test. Only supply this if not all genes were tested and are supplied above in DETest or gene_ids. - :param clean_ref: - Whether or not to only retain gene identifiers in RefSets that occur in + :param clean_ref: Whether or not to only retain gene identifiers in RefSets that occur in the background set of identifiers supplied here through all_ids. - :param capital: - Make all gene IDs captial. + :param capital: Make all gene IDs captial. """ return Enrich( ref=ref, det=det, - pval=pval, + scores=scores, gene_ids=gene_ids, - de_threshold=de_threshold, + threshold=threshold, incl_all_zero=incl_all_zero, all_ids=all_ids, clean_ref=clean_ref, @@ -252,9 +248,9 @@ def __init__( self, ref: RefSets, det: Union[_DifferentialExpressionTest, None], - pval: Union[np.array, None], - gene_ids: Union[list, None], - de_threshold, + scores: Union[np.array, None], + gene_ids: Union[list, np.ndarray, None], + threshold, incl_all_zero, all_ids, clean_ref, @@ -263,6 +259,8 @@ def __init__( self._n_overlaps = None self._pval_enrich = None self._qval_enrich = None + if isinstance(gene_ids, list): + gene_ids = np.asarray(gene_ids) # Load multiple-testing-corrected differential expression # p-values from differential expression output. if det is not None: @@ -273,8 +271,8 @@ def __init__( idx_not_all_zero = np.where(np.logical_not(det.summary()["zero_mean"].values))[0] self._qval_de = det.qval[idx_not_all_zero] self._gene_ids = det.gene_ids[idx_not_all_zero] - elif pval is not None and gene_ids is not None: - self._qval_de = np.asarray(pval) + elif scores is not None and gene_ids is not None: + self._qval_de = np.asarray(scores) self._gene_ids = gene_ids else: raise ValueError('Supply either DETest or pval and gene_ids to Enrich().') @@ -282,7 +280,7 @@ def __init__( # Select significant genes based on user defined threshold. if any([x is np.nan for x in self._gene_ids]): idx_notnan = np.where([x is not np.nan for x in self._gene_ids])[0] - logger.info( + logging.getLogger("diffxpy").info( " Discarded %i nan gene ids, leaving %i genes.", len(self._gene_ids) - len(idx_notnan), len(idx_notnan) @@ -290,7 +288,7 @@ def __init__( self._qval_de = self._qval_de[idx_notnan] self._gene_ids = self._gene_ids[idx_notnan] - self._significant_de = self._qval_de <= de_threshold + self._significant_de = self._qval_de <= threshold self._significant_ids = set(self._gene_ids[np.where(self._significant_de)[0]]) if all_ids is not None: self._all_ids = set(all_ids) @@ -303,7 +301,7 @@ def __init__( self._significant_ids = set([x.upper() for x in self._significant_ids]) # Generate diagnostic statistic of number of possible overlaps in total. - logger.info( + logging.getLogger("diffxpy").info( " %i overlaps found between refset (%i) and provided gene list (%i).", len(set(self._all_ids).intersection(set(ref._genes))), len(ref._genes), @@ -318,7 +316,7 @@ def __init__( # Print if there are empty sets. idx_nonempty = np.where([len(x.genes) > 0 for x in self.RefSets.sets])[0] if len(self.RefSets.sets) - len(idx_nonempty) > 0: - logger.info( + logging.getLogger("diffxpy").info( " Found %i empty sets, removing those.", len(self.RefSets.sets) - len(idx_nonempty) ) @@ -389,7 +387,10 @@ def significant_sets(self, threshold=0.05) -> list: """ Return significant sets from gene set enrichement analysis as an output table. """ - return self.RefSets.subset(idx=np.where(self.qval <= threshold)[0]) + sig_sets = np.where(self.qval <= threshold)[0] + if len(sig_sets) == 0: + logging.getLogger("diffxpy").info("no significant sets found") + return self.RefSets.subset(idx=sig_sets) def significant_set_ids(self, threshold=0.05) -> np.array: """ @@ -424,4 +425,4 @@ def set_summary(self, id: str): :return: Slice of summary table. """ - return self.summary(sort=False).iloc[self.RefSets._ids.index(id), :] + return self.summary(sort=False).iloc[self.RefSets._ids.tolist().index(id), :] diff --git a/diffxpy/testing/__init__.py b/diffxpy/testing/__init__.py index 6836e26..4ab96ba 100644 --- a/diffxpy/testing/__init__.py +++ b/diffxpy/testing/__init__.py @@ -1,3 +1,2 @@ from .tests import lrt, wald, t_test, rank_test, two_sample, pairwise, \ versus_rest, partition, continuous_1d -from .utils import design_matrix, coef_names \ No newline at end of file diff --git a/diffxpy/testing/det.py b/diffxpy/testing/det.py index 3764811..412a15d 100644 --- a/diffxpy/testing/det.py +++ b/diffxpy/testing/det.py @@ -2,6 +2,7 @@ import logging from typing import Union, Dict, Tuple, List, Set import pandas as pd +from random import sample import numpy as np import xarray as xr @@ -14,6 +15,7 @@ except ImportError: anndata = None +from batchglm.xarray_sparse.base import SparseXArrayDataArray from batchglm.models.glm_nb import Model as GeneralizedLinearModel from ..stats import stats @@ -247,7 +249,7 @@ def summary(self, **kwargs) -> pd.DataFrame: def _threshold_summary( self, - res, + res: pd.DataFrame, qval_thres=None, fc_upper_thres=None, fc_lower_thres=None, @@ -255,13 +257,20 @@ def _threshold_summary( ) -> pd.DataFrame: """ Reduce differential expression results into an output table with desired thresholds. + + :param res: Unfiltered summary table. + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Filtered summary table. """ assert fc_lower_thres > 0 if fc_lower_thres is not None else True, "supply positive fc_lower_thres" assert fc_upper_thres > 0 if fc_upper_thres is not None else True, "supply positive fc_upper_thres" if qval_thres is not None: qvals = res['qval'].values - qval_include = np.isnan(qvals) == False + qval_include = np.logical_not(np.isnan(qvals)) qval_include[qval_include] = qvals[qval_include] <= qval_thres res = res.iloc[qval_include, :] @@ -287,12 +296,13 @@ def plot_volcano( alpha=0.05, min_fc=1, size=20, - highlight_ids: List = [], + highlight_ids: Union[List, Tuple] = (), highlight_size: float = 30, highlight_col: str = "red", show: bool = True, save: Union[str, None] = None, - suffix: str = "_volcano.png" + suffix: str = "_volcano.png", + return_axs: bool = False ): """ Returns a volcano plot of p-value vs. log fold change @@ -308,12 +318,13 @@ def plot_volcano( the points below the threshold are colored in grey. :param size: Size of points. :param highlight_ids: Genes to highlight in volcano plot. - :param highlight_ids: Size of points of genes to highlight in volcano plot. - :param highlight_ids: Color of points of genes to highlight in volcano plot. + :param highlight_size: Size of points of genes to highlight in volcano plot. + :param highlight_col: Color of points of genes to highlight in volcano plot. :param show: Whether (if save is not None) and where (save indicates dir and file stem) to display plot. :param save: Path+file name stem to save plots to. File will be save+suffix. Does not save if save is None. :param suffix: Suffix for file name to save plot to. Also use this to set the file type. + :param return_axs: Whether to return axis objects. :return: Tuple of matplotlib (figure, axis) """ @@ -322,7 +333,7 @@ def plot_volcano( plt.ioff() - if corrected_pval == True: + if corrected_pval: neg_log_pvals = - self.log10_qval_clean(log10_threshold=log10_p_threshold) else: neg_log_pvals = - self.log10_pval_clean(log10_threshold=log10_p_threshold) @@ -346,8 +357,8 @@ def plot_volcano( palette={True: "orange", False: "black"}) highlight_ids_found = np.array([x in self.gene_ids for x in highlight_ids]) - highlight_ids_clean = [highlight_ids[i] for i in np.where(highlight_ids_found == True)[0]] - highlight_ids_not_found = [highlight_ids[i] for i in np.where(highlight_ids_found == False)[0]] + highlight_ids_clean = [highlight_ids[i] for i in np.where(highlight_ids_found)[0]] + highlight_ids_not_found = [highlight_ids[i] for i in np.where(np.logical_not(highlight_ids_found))[0]] if len(highlight_ids_not_found) > 0: logger.warning("not all highlight_ids were found in data set: ", ", ".join(highlight_ids_not_found)) @@ -355,8 +366,8 @@ def plot_volcano( neg_log_pvals_highlights = np.zeros([len(highlight_ids_clean)]) logfc_highlights = np.zeros([len(highlight_ids_clean)]) is_highlight = np.zeros([len(highlight_ids_clean)]) - for i,id in enumerate(highlight_ids_clean): - idx = np.where(self.gene_ids == id)[0] + for i, id_i in enumerate(highlight_ids_clean): + idx = np.where(self.gene_ids == id_i)[0] neg_log_pvals_highlights[i] = neg_log_pvals[idx] logfc_highlights[i] = logfc[idx] @@ -365,8 +376,7 @@ def plot_volcano( legend=False, s=highlight_size, palette={0: highlight_col}) - - if corrected_pval == True: + if corrected_pval: ax.set(xlabel="log2FC", ylabel='-log10(corrected p-value)') else: ax.set(xlabel="log2FC", ylabel='-log10(p-value)') @@ -380,7 +390,10 @@ def plot_volcano( plt.close(fig) - return ax + if return_axs: + return ax + else: + return def plot_ma( self, @@ -389,12 +402,13 @@ def plot_ma( min_mean=1e-4, alpha=0.05, size=20, - highlight_ids: List = [], + highlight_ids: Union[List, Tuple] = (), highlight_size: float = 30, highlight_col: str = "red", show: bool = True, save: Union[str, None] = None, - suffix: str = "_my_plot.png" + suffix: str = "_ma_plot.png", + return_axs: bool = False ): """ Returns an MA plot of mean expression vs. log fold change with significance @@ -411,13 +425,13 @@ def plot_ma( non-significant. The corresponding points are colored in grey. :param size: Size of points. :param highlight_ids: Genes to highlight in volcano plot. - :param highlight_ids: Size of points of genes to highlight in volcano plot. - :param highlight_ids: Color of points of genes to highlight in volcano plot. + :param highlight_size: Size of points of genes to highlight in volcano plot. + :param highlight_col: Color of points of genes to highlight in volcano plot. :param show: Whether (if save is not None) and where (save indicates dir and file stem) to display plot. :param save: Path+file name stem to save plots to. File will be save+suffix. Does not save if save is None. :param suffix: Suffix for file name to save plot to. Also use this to set the file type. - + :param return_axs: Whether to return axis objects. :return: Tuple of matplotlib (figure, axis) """ @@ -457,8 +471,8 @@ def plot_ma( palette={True: "orange", False: "black"}) highlight_ids_found = np.array([x in self.gene_ids for x in highlight_ids]) - highlight_ids_clean = [highlight_ids[i] for i in np.where(highlight_ids_found == True)[0]] - highlight_ids_not_found = [highlight_ids[i] for i in np.where(highlight_ids_found == False)[0]] + highlight_ids_clean = [highlight_ids[i] for i in np.where(highlight_ids_found)[0]] + highlight_ids_not_found = [highlight_ids[i] for i in np.where(np.logical_not(highlight_ids_found))[0]] if len(highlight_ids_not_found) > 0: logger.warning("not all highlight_ids were found in data set: ", ", ".join(highlight_ids_not_found)) @@ -466,8 +480,8 @@ def plot_ma( ave_highlights = np.zeros([len(highlight_ids_clean)]) logfc_highlights = np.zeros([len(highlight_ids_clean)]) is_highlight = np.zeros([len(highlight_ids_clean)]) - for i,id in enumerate(highlight_ids_clean): - idx = np.where(self.gene_ids == id)[0] + for i, id_i in enumerate(highlight_ids_clean): + idx = np.where(self.gene_ids == id_i)[0] ave_highlights[i] = ave[idx] logfc_highlights[i] = logfc[idx] @@ -488,7 +502,10 @@ def plot_ma( plt.close(fig) plt.ion() - return ax + if return_axs: + return ax + else: + return class _DifferentialExpressionTestSingle(_DifferentialExpressionTest, metaclass=abc.ABCMeta): @@ -509,6 +526,12 @@ def summary( ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ assert self.gene_ids is not None @@ -517,7 +540,8 @@ def summary( "pval": self.pval, "qval": self.qval, "log2fc": self.log2_fold_change(), - "mean": self.mean + "mean": self.mean, + "zero_mean": self.mean == 0 }) return res @@ -720,11 +744,22 @@ def scales(self): return retval - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ res = super().summary(**kwargs) res["grad"] = self.full_model_gradient.data @@ -747,6 +782,7 @@ class DifferentialExpressionTestWald(_DifferentialExpressionTestSingle): """ model_estim: _Estimation + sample_description: pd.DataFrame coef_loc_totest: np.ndarray theta_mle: np.ndarray theta_sd: np.ndarray @@ -756,16 +792,21 @@ class DifferentialExpressionTestWald(_DifferentialExpressionTestSingle): def __init__( self, model_estim: _Estimation, - col_indices: np.ndarray + col_indices: np.ndarray, + noise_model: str, + sample_description: pd.DataFrame ): """ :param model_estim: - :param cold_index: indices of coefs to test + :param col_indices: indices of coefs to test """ super().__init__() + self.sample_description = sample_description self.model_estim = model_estim self.coef_loc_totest = col_indices + self.noise_model = noise_model + self._store_ols = None try: if model_estim._error_codes is not None: @@ -859,11 +900,22 @@ def _test(self): theta0=0 ) - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ res = super().summary(**kwargs) res["grad"] = self.model_gradient.data @@ -889,14 +941,29 @@ def summary(self, qval_thres=None, fc_upper_thres=None, return res - def plot_vs_ttest(self, log10=False): + def plot_vs_ttest( + self, + log10=False, + return_axs: bool = False + ): + """ + Normalizes data by size factors if any were used in model. + + :param log10: + :param return_axs: Whether to return axis objects. + + :return: + """ 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]) + # Normalize by size factors that were used in regression. + sf = np.broadcast_to(np.expand_dims(self.model_estim.size_factors, axis=1), + shape=self.model_estim.X.shape) ttest = t_test( - data=self.model_estim.X, + data=self.model_estim.X.multiply(1 / sf, copy=True), grouping=grouping, gene_names=self.gene_ids, ) @@ -913,7 +980,622 @@ def plot_vs_ttest(self, log10=False): ax.set(xlabel="t-test", ylabel='wald test') - return fig, ax + if return_axs: + return ax + else: + return + + def plot_comparison_ols_coef( + self, + size=20, + show: bool = True, + save: Union[str, None] = None, + suffix: str = "_ols_comparison_coef.png", + ncols=3, + row_gap=0.3, + col_gap=0.25, + return_axs: bool = False + ): + """ + Plot location model coefficients of inferred model against those obtained from an OLS model. + + Red line shown is the identity line. + Note that this comparison only seems to be useful if the covariates are zero centred. This is + especially important for continuous covariates. + + :param size: Size of points. + :param show: Whether (if save is not None) and where (save indicates dir and file stem) to display plot. + :param save: Path+file name stem to save plots to. + File will be save+suffix. Does not save if save is None. + :param suffix: Suffix for file name to save plot to. Also use this to set the file type. + :param ncols: Number of columns in plot grid if multiple genes are plotted. + :param row_gap: Vertical gap between panel rows relative to panel height. + :param col_gap: Horizontal gap between panel columns relative to panel width. + :param return_axs: Whether to return axis objects. + + :return: Matplotlib axis objects. + """ + import seaborn as sns + import matplotlib.pyplot as plt + from matplotlib import gridspec + from matplotlib import rcParams + from batchglm.api.models.glm_norm import Estimator, InputData + + # Run OLS model fit to have comparison coefficients. + if self._store_ols is None: + input_data_ols = InputData.new( + data=self.model_estim.input_data.data, + design_loc=self.model_estim.input_data.design_loc, + design_scale=self.model_estim.input_data.design_scale[:, [0]], + constraints_loc=self.model_estim.input_data.constraints_loc, + constraints_scale=self.model_estim.input_data.constraints_scale[[0], [0]], + size_factors=self.model_estim.input_data.size_factors, + feature_names=self.model_estim.input_data.features, + ) + estim_ols = Estimator( + input_data=input_data_ols, + init_model=None, + init_a="standard", + init_b="standard", + dtype=self.model_estim.a_var.dtype + ) + estim_ols.initialize() + store_ols = estim_ols.finalize() + self._store_ols = store_ols + else: + store_ols = self._store_ols + + # Prepare parameter summary of both model fits. + par_loc = self.model_estim.input_data.data.coords["design_loc_params"].values + + a_var_ols = store_ols.a_var.values + a_var_ols[1:, :] = (a_var_ols[1:, :] + a_var_ols[[0], :]) / a_var_ols[[0], :] + + a_var_user = self.model_estim.a_var.values + # Translate coefficients from both fits to be multiplicative in identity space. + if self.noise_model == "nb": + a_var_user = np.exp(a_var_user) # self.model_estim.inverse_link_loc(a_var_user) + elif self.noise_model == "norm": + a_var_user[1:, :] = (a_var_user[1:, :] + a_var_user[[0], :]) / a_var_user[[0], :] + else: + raise ValueError("noise model %s not yet supported for plot_comparison_ols" % self.noise_model) + + summaries_fits = [ + pd.DataFrame({ + "user": a_var_user[i, :], + "ols": a_var_ols[i, :], + "coef": par_loc[i] + }) for i in range(self.model_estim.a_var.shape[0]) + ] + + plt.ioff() + nrows = len(par_loc) // ncols + int((len(par_loc) % ncols) > 0) + + gs = gridspec.GridSpec( + nrows=nrows, + ncols=ncols, + hspace=row_gap, + wspace=col_gap + ) + fig = plt.figure( + figsize=( + ncols * rcParams['figure.figsize'][0], # width in inches + nrows * rcParams['figure.figsize'][1] * (1 + row_gap) # height in inches + ) + ) + + axs = [] + for i, par_i in enumerate(par_loc): + ax = plt.subplot(gs[i]) + axs.append(ax) + + x = summaries_fits[i]["user"].values + y = summaries_fits[i]["ols"].values + + sns.scatterplot( + x=x, + y=y, + ax=ax, + s=size + ) + sns.lineplot( + x=np.array([np.min([np.min(x), np.min(y)]), np.max([np.max(x), np.max(y)])]), + y=np.array([np.min([np.min(x), np.min(y)]), np.max([np.max(x), np.max(y)])]), + ax=ax, + color="red", + legend=False + ) + ax.set(xlabel="user supplied model", ylabel="OLS model") + title_i = par_loc[i] + " (R=" + str(np.round(np.corrcoef(x, y)[0, 1], 3)) + ")" + ax.set_title(title_i) + + # Save, show and return figure. + if save is not None: + plt.savefig(save + suffix) + + if show: + plt.show() + + plt.close(fig) + plt.ion() + + if return_axs: + return axs + else: + return + + def plot_comparison_ols_pred( + self, + size=20, + log1p_transform: bool = True, + show: bool = True, + save: Union[str, None] = None, + suffix: str = "_ols_comparison_pred.png", + row_gap=0.3, + col_gap=0.25, + return_axs: bool = False + ): + """ + Compare location model prediction of inferred model with one obtained from an OLS model. + + Red line shown is the identity line. + + :param size: Size of points. + :param log1p_transform: Whether to log1p transform the data. + :param show: Whether (if save is not None) and where (save indicates dir and file stem) to display plot. + :param save: Path+file name stem to save plots to. + File will be save+suffix. Does not save if save is None. + :param suffix: Suffix for file name to save plot to. Also use this to set the file type. + :param row_gap: Vertical gap between panel rows relative to panel height. + :param col_gap: Horizontal gap between panel columns relative to panel width. + :param return_axs: Whether to return axis objects. + + :return: Matplotlib axis objects. + """ + import seaborn as sns + import matplotlib.pyplot as plt + from matplotlib import gridspec + from matplotlib import rcParams + from batchglm.api.models.glm_norm import Estimator, InputData + + # Run OLS model fit to have comparison coefficients. + if self._store_ols is None: + input_data_ols = InputData.new( + data=self.model_estim.input_data.data, + design_loc=self.model_estim.input_data.design_loc, + design_scale=self.model_estim.input_data.design_scale[:, [0]], + constraints_loc=self.model_estim.input_data.constraints_loc, + constraints_scale=self.model_estim.input_data.constraints_scale[[0], [0]], + size_factors=self.model_estim.input_data.size_factors, + feature_names=self.model_estim.input_data.features, + ) + estim_ols = Estimator( + input_data=input_data_ols, + init_model=None, + init_a="standard", + init_b="standard", + dtype=self.model_estim.a_var.dtype + ) + estim_ols.initialize() + store_ols = estim_ols.finalize() + self._store_ols = store_ols + else: + store_ols = self._store_ols + + # Prepare parameter summary of both model fits. + plt.ioff() + nrows = 1 + ncols = 2 + + axs = [] + gs = gridspec.GridSpec( + nrows=nrows, + ncols=ncols, + hspace=row_gap, + wspace=col_gap + ) + fig = plt.figure( + figsize=( + ncols * rcParams['figure.figsize'][0], # width in inches + nrows * rcParams['figure.figsize'][1] * (1 + row_gap) # height in inches + ) + ) + + pred_n_cells = sample( + population=list(np.arange(0, self.model_estim.X.shape[0])), + k=np.min([20, self.model_estim.design_loc.shape[0]]) + ) + + if isinstance(self.model_estim.X, SparseXArrayDataArray): + x = np.asarray(self.model_estim.X.X[pred_n_cells, :].todense()).flatten() + else: + x = np.asarray(self.model_estim.X[pred_n_cells, :]).flatten() + + y_user = self.model_estim.inverse_link_loc( + np.matmul(self.model_estim.design_loc[pred_n_cells, :].values, self.model_estim.a_var.values).flatten() + ) + y_ols = store_ols.inverse_link_loc( + np.matmul(store_ols.design_loc[pred_n_cells, :].values, store_ols.a_var.values).flatten() + ) + if log1p_transform: + x = np.log(x+1) + y_user = np.log(y_user + 1) + y_ols = np.log(y_ols + 1) + + y = np.concatenate([y_user, y_ols]) + + summary0_fit = pd.concat([ + pd.DataFrame({ + "observed": y_user, + "predicted": x, + "model": ["user" for i in x] + }), + pd.DataFrame({ + "observed": y_ols, + "predicted": x, + "model": ["OLS" for i in x] + }) + ]) + + ax0 = plt.subplot(gs[0]) + axs.append(ax0) + sns.scatterplot( + x="observed", + y="predicted", + hue="model", + data=summary0_fit, + ax=ax0, + s=size + ) + sns.lineplot( + x=np.array([np.min([np.min(x), np.min(y)]), np.max([np.max(x), np.max(y)])]), + y=np.array([np.min([np.min(x), np.min(y)]), np.max([np.max(x), np.max(y)])]), + ax=ax0, + color="red", + legend=False + ) + ax0.set(xlabel="observed value", ylabel="model") + + summary1_fit = pd.concat([ + pd.DataFrame({ + "dev": y_user-x, + "model": ["user" for i in x] + }), + pd.DataFrame({ + "dev": y_ols-x, + "model": ["OLS" for i in x] + }) + ]) + + ax1 = plt.subplot(gs[1]) + axs.append(ax0) + sns.boxplot( + x="model", + y="dev", + data=summary1_fit, + ax=ax1 + ) + ax1.set(xlabel="model", ylabel="deviation from observations") + + # Save, show and return figure. + if save is not None: + plt.savefig(save + suffix) + + if show: + plt.show() + + plt.close(fig) + plt.ion() + + if return_axs: + return axs + else: + return + + def _assemble_gene_fits( + self, + gene_names: Tuple, + covariate_x: str, + covariate_hue: Union[None, str], + log1p_transform: bool, + incl_fits: bool + ): + """ + Prepare data for gene-wise model plots. + + :param gene_names: Genes to generate plots for. + :param covariate_x: Covariate in location model to partition x-axis by. + :param covariate_hue: Covariate in location model to stack boxplots by. + :param log1p_transform: Whether to log transform observations + before estimating the distribution with boxplot. Model estimates are adjusted accordingly. + :param incl_fits: Whether to include fits in plot. + :return summaries_genes: List with data frame for seabron in it. + """ + + summaries_genes = [] + for i, g in enumerate(gene_names): + assert g in self.model_estim.features, "gene %g not found" % g + g_idx = self.model_estim.features.tolist().index(g) + # Raw data for boxplot: + y = self.model_estim.X[:, g_idx] + # Model fits: + loc = self.model_estim.location[:, g_idx] + scale = self.model_estim.scale[:, g_idx] + if self.noise_model == "nb": + yhat = np.random.negative_binomial( + n=scale, + p=1 - loc / (scale + loc) + ) + elif self.noise_model == "norm": + yhat = np.random.normal( + loc=loc, + scale=scale + ) + else: + raise ValueError("noise model %s not yet supported for plot_gene_fits" % self.noise_model) + + # Transform observed data: + if log1p_transform: + y = np.log(y + 1) + yhat = np.log(yhat + 1) + + # Build DataFrame which contains all information for raw data: + summary_raw = pd.DataFrame({"y": y, "data": "obs"}) + if incl_fits: + summary_fit = pd.DataFrame({"y": yhat, "data": "fit"}) + if covariate_x is not None: + assert self.sample_description is not None, "sample_description was not provided to test.wald()" + if covariate_x in self.sample_description.columns: + summary_raw["x"] = self.sample_description[covariate_x].values.astype(str) + if incl_fits: + summary_fit["x"] = self.sample_description[covariate_x].values.astype(str) + else: + raise ValueError("covariate_x=%s not found in location model" % covariate_x) + else: + summary_raw["x"] = " " + if incl_fits: + summary_fit["x"] = " " + + if covariate_hue is not None: + assert self.sample_description is not None, "sample_description was not provided to test.wald()" + if covariate_hue in self.sample_description.columns: + if incl_fits: + summary_raw["hue"] = [str(x)+"_obs" for x in self.sample_description[covariate_hue].values] + summary_fit["hue"] = [str(x)+"_fit" for x in self.sample_description[covariate_hue].values] + else: + summary_raw["hue"] = self.sample_description[covariate_hue].values + else: + raise ValueError("covariate_x=%s not found in location model" % covariate_x) + else: + summary_raw["hue"] = "obs" + if incl_fits: + summary_fit["hue"] = "fit" + + if incl_fits: + summaries = pd.concat([summary_raw, summary_fit]) + else: + summaries = summary_raw + summaries.x = pd.Categorical(summaries.x, ordered=True) + summaries.hue = pd.Categorical(summaries.hue, ordered=True) + + summaries_genes.append(summaries) + + return summaries_genes + + def plot_gene_fits_boxplots( + self, + gene_names: Tuple, + covariate_x: str = None, + covariate_hue: str = None, + log1p_transform: bool = False, + incl_fits: bool = True, + show: bool = True, + save: Union[str, None] = None, + suffix: str = "_genes_boxplot.png", + ncols=3, + row_gap=0.3, + col_gap=0.25, + xtick_rotation=0, + legend: bool = True, + return_axs: bool = False, + **kwargs + ): + """ + Plot gene-wise model fits and observed distribution by covariates. + + Use this to inspect fitting performance on individual genes. + + :param gene_names: Genes to generate plots for. + :param covariate_x: Covariate in location model to partition x-axis by. + :param covariate_hue: Covariate in location model to stack boxplots by. + :param log1p_transform: Whether to log transform observations + before estimating the distribution with boxplot. Model estimates are adjusted accordingly. + :param incl_fits: Whether to include fits in plot. + :param show: Whether (if save is not None) and where (save indicates dir and file stem) to display plot. + :param save: Path+file name stem to save plots to. + File will be save+suffix. Does not save if save is None. + :param suffix: Suffix for file name to save plot to. Also use this to set the file type. + :param ncols: Number of columns in plot grid if multiple genes are plotted. + :param row_gap: Vertical gap between panel rows relative to panel height. + :param col_gap: Horizontal gap between panel columns relative to panel width. + :param xtick_rotation: Angle to rotate x-ticks by. + :param legend: Whether to show legend. + :param return_axs: Whether to return axis objects. + + :return: Matplotlib axis objects. + """ + import seaborn as sns + import matplotlib.pyplot as plt + from matplotlib import gridspec + from matplotlib import rcParams + + plt.ioff() + nrows = len(gene_names) // ncols + int((len(gene_names) % ncols) > 0) + + gs = gridspec.GridSpec( + nrows=nrows, + ncols=ncols, + hspace=row_gap, + wspace=col_gap + ) + fig = plt.figure( + figsize=( + ncols * rcParams['figure.figsize'][0], # width in inches + nrows * rcParams['figure.figsize'][1] * (1 + row_gap) # height in inches + ) + ) + + axs = [] + summaries = self._assemble_gene_fits( + gene_names=gene_names, + covariate_x=covariate_x, + covariate_hue=covariate_hue, + log1p_transform=log1p_transform, + incl_fits=incl_fits + ) + for i, g in enumerate(gene_names): + ax = plt.subplot(gs[i]) + axs.append(ax) + + if log1p_transform: + ylabel = "log1p expression" + else: + ylabel = "expression" + + sns.boxplot( + x="x", + y="y", + hue="hue", + data=summaries[i], + ax=ax, + **kwargs + ) + + ax.set(xlabel="covariate", ylabel=ylabel) + ax.set_title(g) + if not legend: + ax.legend_.remove() + + plt.xticks(rotation=xtick_rotation) + + # Save, show and return figure. + if save is not None: + plt.savefig(save + suffix) + + if show: + plt.show() + + plt.close(fig) + plt.ion() + + if return_axs: + return axs + else: + return + + def plot_gene_fits_violins( + self, + gene_names: Tuple, + covariate_x: str = None, + log1p_transform: bool = False, + show: bool = True, + save: Union[str, None] = None, + suffix: str = "_genes_violin.png", + ncols=3, + row_gap=0.3, + col_gap=0.25, + xtick_rotation=0, + return_axs: bool = False, + **kwargs + ): + """ + Plot gene-wise model fits and observed distribution by covariates as violins. + + Use this to inspect fitting performance on individual genes. + + :param gene_names: Genes to generate plots for. + :param covariate_x: Covariate in location model to partition x-axis by. + :param log1p_transform: Whether to log transform observations + before estimating the distribution with boxplot. Model estimates are adjusted accordingly. + :param show: Whether (if save is not None) and where (save indicates dir and file stem) to display plot. + :param save: Path+file name stem to save plots to. + File will be save+suffix. Does not save if save is None. + :param suffix: Suffix for file name to save plot to. Also use this to set the file type. + :param ncols: Number of columns in plot grid if multiple genes are plotted. + :param row_gap: Vertical gap between panel rows relative to panel height. + :param col_gap: Horizontal gap between panel columns relative to panel width. + :param xtick_rotation: Angle to rotate x-ticks by. + :param return_axs: Whether to return axis objects. + + :return: Matplotlib axis objects. + """ + import seaborn as sns + import matplotlib.pyplot as plt + from matplotlib import gridspec + from matplotlib import rcParams + + plt.ioff() + nrows = len(gene_names) // ncols + int((len(gene_names) % ncols) > 0) + + gs = gridspec.GridSpec( + nrows=nrows, + ncols=ncols, + hspace=row_gap, + wspace=col_gap + ) + fig = plt.figure( + figsize=( + ncols * rcParams['figure.figsize'][0], # width in inches + nrows * rcParams['figure.figsize'][1] * (1 + row_gap) # height in inches + ) + ) + + axs = [] + summaries = self._assemble_gene_fits( + gene_names=gene_names, + covariate_x=covariate_x, + covariate_hue=None, + log1p_transform=log1p_transform, + incl_fits=True + ) + for i, g in enumerate(gene_names): + ax = plt.subplot(gs[i]) + axs.append(ax) + + if log1p_transform: + ylabel = "log1p expression" + else: + ylabel = "expression" + + sns.violinplot( + x="x", + y="y", + hue="data", + split=True, + data=summaries[i], + ax=ax, + **kwargs + ) + + ax.set(xlabel="covariate", ylabel=ylabel) + ax.set_title(g) + + plt.xticks(rotation=xtick_rotation) + + # Save, show and return figure. + if save is not None: + plt.savefig(save + suffix) + + if show: + plt.show() + + plt.close(fig) + plt.ion() + + if return_axs: + return axs + else: + return class DifferentialExpressionTestTT(_DifferentialExpressionTestSingle): @@ -921,9 +1603,18 @@ class DifferentialExpressionTestTT(_DifferentialExpressionTestSingle): Single t-test test per gene. """ - def __init__(self, data, grouping, gene_names, is_logged): + def __init__( + self, + data, + sample_description: pd.DataFrame, + grouping, + gene_names, + is_logged, + is_sig_zerovar: bool = True + ): super().__init__() self._X = data + self.sample_description = sample_description self.grouping = grouping self._gene_names = np.asarray(gene_names) @@ -957,6 +1648,16 @@ def __init__(self, data, grouping, gene_names, is_logged): n0=x0.shape[0], n1=x1.shape[0] ) + pval[np.where(np.logical_and( + np.logical_and(mean_x0 == mean_x1, self._mean > 0), + np.logical_not(self._var_geq_zero) + ))[0]] = 1.0 + if is_sig_zerovar: + pval[np.where(np.logical_and( + mean_x0 != mean_x1, + np.logical_not(self._var_geq_zero) + ))[0]] = 0.0 + self._pval = pval if is_logged: @@ -969,20 +1670,20 @@ def __init__(self, data, grouping, gene_names, is_logged): # This is the default which can be changed and can be changed # via DIFFXPY_TREAT_ZEROVAR_TT_AS_SIG. pval[np.where(np.logical_and(np.logical_and( - self._var_geq_zero == False, - self._ave_nonzero == True), + np.logical_not(self._var_geq_zero), + self._ave_nonzero), np.abs(self._logfc) < np.nextafter(0, 1) ))] = 0 if pkg_constants.DE_TREAT_ZEROVAR_TT_AS_SIG: pval[np.where(np.logical_and(np.logical_and( - self._var_geq_zero == False, - self._ave_nonzero == True), + np.logical_not(self._var_geq_zero), + self._ave_nonzero), np.abs(self._logfc) >= np.nextafter(0, 1) ))] = 1 else: pval[np.where(np.logical_and(np.logical_and( - self._var_geq_zero == False, - self._ave_nonzero == True), + np.logical_not(self._var_geq_zero), + self._ave_nonzero), np.abs(self._logfc) >= np.nextafter(0, 1) ))] = 0 @@ -1000,15 +1701,25 @@ def log_fold_change(self, base=np.e, **kwargs): """ return self._logfc / np.log(base) - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ res = super().summary(**kwargs) - res["zero_mean"] = self._ave_nonzero == False - res["zero_variance"] = self._var_geq_zero == False + res["zero_variance"] = np.logical_not(self._var_geq_zero) res = self._threshold_summary( res=res, @@ -1026,9 +1737,18 @@ class DifferentialExpressionTestRank(_DifferentialExpressionTestSingle): Single rank test per gene (Mann-Whitney U test). """ - def __init__(self, data, grouping, gene_names, is_logged): + def __init__( + self, + data, + sample_description: pd.DataFrame, + grouping, + gene_names, + is_logged, + is_sig_zerovar: bool = True + ): super().__init__() self._X = data + self.sample_description = sample_description self.grouping = grouping self._gene_names = np.asarray(gene_names) @@ -1064,6 +1784,15 @@ def __init__(self, data, grouping, gene_names, is_logged): x0=x0.X[:, idx_run].toarray(), x1=x1.X[:, idx_run].toarray() ) + pval[np.where(np.logical_and( + np.logical_and(mean_x0 == mean_x1, self._mean > 0), + np.logical_not(self._var_geq_zero) + ))[0]] = 1.0 + if is_sig_zerovar: + pval[np.where(np.logical_and( + mean_x0 != mean_x1, + np.logical_not(self._var_geq_zero) + ))[0]] = 0.0 self._pval = pval @@ -1089,13 +1818,25 @@ def log_fold_change(self, base=np.e, **kwargs): else: return self._logfc / np.log(base) - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ res = super().summary(**kwargs) + res["zero_variance"] = np.logical_not(self._var_geq_zero) res = self._threshold_summary( res=res, @@ -1160,7 +1901,7 @@ def _correction(self, method): return qvals elif self._correction_type.lower() == "by_test": qvals = np.apply_along_axis( - func1d=lambda pvals: correction.correct(pvals=pvals, method=method), + func1d=lambda pv: correction.correct(pvals=pv, method=method), axis=-1, arr=self.pval, ) @@ -1219,7 +1960,7 @@ def __init__(self, gene_ids, pval, logfc, ave, groups, tests, correction_type: s self.groups = list(np.asarray(groups)) self._tests = tests - q = self.qval + _ = self.qval @property def gene_ids(self) -> np.ndarray: @@ -1318,12 +2059,18 @@ def log10_qval_pair_clean(self, group1, group2, log10_threshold=-30): log10_qval_clean = np.clip(log10_qval_clean, log10_threshold, 0, log10_qval_clean) return log10_qval_clean - def log_fold_change_pair(self, group1, group2, base=np.e): + def log_fold_change_pair( + self, + group1, + group2, + base=np.e + ): """ Get log fold changes of the comparison of group1 and group2. :param group1: Identifier of first group of observations in pair-wise comparison. :param group2: Identifier of second group of observations in pair-wise comparison. + :param base: Base of logarithm. :return: log fold changes """ assert self._logfc is not None @@ -1331,11 +2078,22 @@ def log_fold_change_pair(self, group1, group2, base=np.e): self._check_groups(group1, group2) return self.log_fold_change(base=base)[self.groups.index(group1), self.groups.index(group2), :] - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ res = super().summary(**kwargs) @@ -1349,15 +2107,24 @@ def summary(self, qval_thres=None, fc_upper_thres=None, return res - def summary_pair(self, group1, group2, - qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary_pair( + self, + group1, + group2, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. :param group1: Identifier of first group of observations in pair-wise comparison. :param group2: Identifier of second group of observations in pair-wise comparison. + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. :return: pandas.DataFrame with the following columns: - gene: the gene id's @@ -1410,8 +2177,8 @@ def __init__(self, model_estim: _Estimation, grouping, groups, correction_type: self._logfc = None # Call tests in constructor. - p = self.pval - q = self.qval + _ = self.pval + _ = self.qval def _test(self, **kwargs): groups = self.groups @@ -1503,11 +2270,22 @@ def log_fold_change_pair(self, group1, group2, base=np.e): self._check_groups(group1, group2) return self.log_fold_change(base=base)[self.groups.index(group1), self.groups.index(group2), :] - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ res = super().summary(**kwargs) @@ -1521,13 +2299,24 @@ def summary(self, qval_thres=None, fc_upper_thres=None, return res - def summary_pair(self, group1, group2, - qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary_pair( + self, + group1, + group2, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + :param group1: Identifier of first group of observations in pair-wise comparison. + :param group2: Identifier of second group of observations in pair-wise comparison. + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. :return: pandas.DataFrame with the following columns: - gene: the gene id's @@ -1556,6 +2345,7 @@ def summary_pair(self, group1, group2, return res + class DifferentialExpressionTestZTestLazy(_DifferentialExpressionTestMulti): """ Pairwise unit_test between more than 2 groups per gene with lazy evaluation. @@ -1619,13 +2409,13 @@ def _test(self, **kwargs): """ pass - def _test_pairs(self, groups0, groups1, **kwargs): + def _test_pairs(self, groups0, groups1): num_features = self.model_estim.X.shape[1] pvals = np.tile(np.NaN, [len(groups0), len(groups1), num_features]) - for i,g0 in enumerate(groups0): - for j,g1 in enumerate(groups1): + for i, g0 in enumerate(groups0): + for j, g1 in enumerate(groups1): if g0 != g1: pvals[i, j] = stats.two_coef_z_test( theta_mle0=self._theta_mle[g0], @@ -1686,12 +2476,25 @@ def log_fold_change(self, base=np.e, **kwargs): """ pass - def summary(self, qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary( + self, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ + Summarize differential expression results into an output table. + This function is not available in lazy results evaluation as it would require all pairwise tests to be performed. + + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. + :return: Summary table of differential expression test. """ pass @@ -1708,7 +2511,7 @@ def _check_groups(self, groups0, groups1): raise ValueError('groups1 element '+str(g)+' not recognized') def _groups_idx(self, groups): - if isinstance(groups, list)==False: + if not isinstance(groups, list): groups = [groups] return np.array([self.groups.index(x) for x in groups]) @@ -1782,25 +2585,35 @@ def log_fold_change_pairs(self, groups0=None, groups1=None, base=np.e): num_features = self._theta_mle.shape[1] logfc = np.zeros(shape=(len(groups0), len(groups1), num_features)) - for i,g0 in enumerate(groups0): - for j,g1 in enumerate(groups1): - logfc[i,j,:] = self._theta_mle[g0,:].values - self._theta_mle[g1,:].values + for i, g0 in enumerate(groups0): + for j, g1 in enumerate(groups1): + logfc[i, j, :] = self._theta_mle[g0, :].values - self._theta_mle[g1, :].values if base == np.e: return logfc else: return logfc / np.log(base) - def summary_pair(self, group0, group1, - qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary_pair( + self, + group0, + group1, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results of single pairwose comparison into an output table. - :param group0: Firt group in pair-wise comparison. - :param group1: Second group in pair-wise comparison. + :param group0: First set of groups in pair-wise comparison. + :param group1: Second set of groups in pair-wise comparison. + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. :return: pandas.DataFrame with the following columns: - gene: the gene id's @@ -1836,16 +2649,26 @@ def summary_pair(self, group0, group1, return res - def summary_pairs(self, groups0, groups1=None, - qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary_pairs( + self, + groups0, + groups1=None, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None, + **kwargs + ) -> pd.DataFrame: """ Summarize differential expression results of a set of pairwise comparisons into an output table. :param groups0: First set of groups in pair-wise comparison. :param groups1: Second set of groups in pair-wise comparison. + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. :return: pandas.DataFrame with the following columns: - gene: the gene id's @@ -1873,8 +2696,8 @@ def summary_pairs(self, groups0, groups1=None, res = pd.DataFrame({ "gene": self.gene_ids, - "pval": np.min(pval, axis=(0,1)), - "qval": np.min(qval, axis=(0,1)), + "pval": np.min(pval, axis=(0, 1)), + "qval": np.min(qval, axis=(0, 1)), "log2fc": np.asarray(logfc), "mean": np.asarray(self.mean) }) @@ -1904,7 +2727,7 @@ def __init__(self, gene_ids, pval, logfc, ave, groups, tests, correction_type: s self.groups = list(np.asarray(groups)) self._tests = tests - q = self.qval + _ = self.qval @property def tests(self): @@ -1921,7 +2744,7 @@ def gene_ids(self) -> np.ndarray: return self._gene_ids @property - def X(self) -> np.ndarray: + def X(self) -> Union[np.ndarray, None]: return None def log_fold_change(self, base=np.e, **kwargs): @@ -1964,13 +2787,22 @@ def summary(self, qval_thres=None, fc_upper_thres=None, return res - def summary_group(self, group, - qval_thres=None, fc_upper_thres=None, - fc_lower_thres=None, mean_thres=None, - **kwargs) -> pd.DataFrame: + def summary_group( + self, + group, + qval_thres=None, + fc_upper_thres=None, + fc_lower_thres=None, + mean_thres=None + ) -> pd.DataFrame: """ Summarize differential expression results into an output table. + :param group: + :param qval_thres: Upper bound of corrected p-values for gene to be included. + :param fc_upper_thres: Upper bound of fold-change for gene to be included. + :param fc_lower_thres: Lower bound of fold-change p-values for gene to be included. + :param mean_thres: Lower bound of average expression for gene to be included. :return: pandas.DataFrame with the following columns: - gene: the gene id's @@ -2014,7 +2846,7 @@ def __init__(self, partitions, tests, ave, correction_type: str = "by_test"): self._logfc = np.expand_dims(np.vstack([x.log_fold_change() for x in tests]), axis=0) self._mean = ave - q = self.qval + _ = self.qval @property def gene_ids(self) -> np.ndarray: @@ -2081,14 +2913,16 @@ def __init__( de_test: _DifferentialExpressionTestSingle, model_estim: _Estimation, size_factors: np.ndarray, - continuous_coords: str, - spline_coefs: list + continuous_coords: np.ndarray, + spline_coefs: list, + noise_model: str ): self._de_test = de_test self._model_estim = model_estim self._size_factors = size_factors self._continuous_coords = continuous_coords self._spline_coefs = spline_coefs + self.noise_model = noise_model @property def gene_ids(self) -> np.ndarray: @@ -2157,8 +2991,7 @@ def log_fold_change(self, base=np.e, genes=None, nonnumeric=False): else: genes = self._idx_genes(genes) - fc = self.max(genes=genes, nonnumeric=nonnumeric) - \ - self.min(genes=genes, nonnumeric=nonnumeric) + fc = self.max(genes=genes, non_numeric=nonnumeric) - self.min(genes=genes, non_numeric=nonnumeric) fc = np.nextafter(0, 1, out=fc, where=fc == 0) return np.log(fc) / np.log(base) @@ -2171,7 +3004,7 @@ def _filter_genes_str(self, genes: list): :return: Filtered list of genes """ genes_found = np.array([x in self.gene_ids for x in genes]) - if any(genes_found == False): + if any(np.logical_not(genes_found)): logger.info("did not find some genes, omitting") genes = genes[genes_found] return genes @@ -2219,75 +3052,75 @@ def _spline_par_loc_idx(self, intercept=True): idx = np.concatenate([np.where([[x == 'Intercept' for x in par_loc_names]])[0], idx]) return idx - def _continuous_model(self, idx, nonnumeric=False): + def _continuous_model(self, idx, non_numeric=False): """ Recover continuous fit for a gene. :param idx: Index of genes to recover fit for. - :param nonnumeric: Whether to include non-numeric covariates in fit. + :param non_numeric: Whether to include non-numeric covariates in fit. :return: Continuuos fit for each cell for given gene. """ idx = np.asarray(idx) - if nonnumeric: + if non_numeric: mu = np.matmul(self._model_estim.design_loc.values, - self._model_estim.par_link_loc[:,idx]) + self._model_estim.par_link_loc[:, idx]) if self._size_factors is not None: mu = mu + self._size_factors else: idx_basis = self._spline_par_loc_idx(intercept=True) - mu = np.matmul(self._model_estim.design_loc[:,idx_basis].values, + mu = np.matmul(self._model_estim.design_loc[:, idx_basis].values, self._model_estim.par_link_loc[idx_basis, idx]) mu = np.exp(mu) return mu - def max(self, genes, nonnumeric=False): + def max(self, genes, non_numeric=False): """ Return maximum fitted expression value by gene. :param genes: Genes for which to return maximum fitted value. - :param nonnumeric: Whether to include non-numeric covariates in fit. + :param non_numeric: Whether to include non-numeric covariates in fit. :return: Maximum fitted expression value by gene. """ genes = self._idx_genes(genes) - return np.array([np.max(self._continuous_model(idx=i, nonnumeric=nonnumeric)) + return np.array([np.max(self._continuous_model(idx=i, non_numeric=non_numeric)) for i in genes]) - def min(self, genes, nonnumeric=False): + def min(self, genes, non_numeric=False): """ Return minimum fitted expression value by gene. :param genes: Genes for which to return maximum fitted value. - :param nonnumeric: Whether to include non-numeric covariates in fit. + :param non_numeric: Whether to include non-numeric covariates in fit. :return: Maximum fitted expression value by gene. """ genes = self._idx_genes(genes) - return np.array([np.min(self._continuous_model(idx=i, nonnumeric=nonnumeric)) + return np.array([np.min(self._continuous_model(idx=i, non_numeric=non_numeric)) for i in genes]) - def argmax(self, genes, nonnumeric=False): + def argmax(self, genes, non_numeric=False): """ Return maximum fitted expression value by gene. :param genes: Genes for which to return maximum fitted value. - :param nonnumeric: Whether to include non-numeric covariates in fit. + :param non_numeric: Whether to include non-numeric covariates in fit. :return: Maximum fitted expression value by gene. """ genes = self._idx_genes(genes) - idx = np.array([np.argmax(self._continuous_model(idx=i, nonnumeric=nonnumeric)) + idx = np.array([np.argmax(self._continuous_model(idx=i, non_numeric=non_numeric)) for i in genes]) return self._continuous_coords[idx] - def argmin(self, genes, nonnumeric=False): + def argmin(self, genes, non_numeric=False): """ Return minimum fitted expression value by gene. :param genes: Genes for which to return maximum fitted value. - :param nonnumeric: Whether to include non-numeric covariates in fit. + :param non_numeric: Whether to include non-numeric covariates in fit. :return: Maximum fitted expression value by gene. """ genes = self._idx_genes(genes) - idx = np.array([np.argmin(self._continuous_model(idx=i, nonnumeric=nonnumeric)) + idx = np.array([np.argmin(self._continuous_model(idx=i, non_numeric=non_numeric)) for i in genes]) return self._continuous_coords[idx] @@ -2297,7 +3130,7 @@ def plot_genes( hue=None, size=1, log=True, - nonnumeric=False, + non_numeric=False, save=None, show=True, ncols=2, @@ -2311,7 +3144,8 @@ def plot_genes( :param genes: Gene IDs to plot. :param hue: Confounder to include in plot. :param size: Point size. - :param nonnumeric: + :param log: Whether to log values. + :param non_numeric: :param save: Path+file name stem to save plots to. File will be save+"_genes.png". Does not save if save is None. :param show: Whether to display plot. @@ -2356,7 +3190,7 @@ def plot_genes( axs.append(ax) y = self.X[:, g] - yhat = self._continuous_model(idx=g, nonnumeric=nonnumeric) + yhat = self._continuous_model(idx=g, non_numeric=non_numeric) if log: y = np.log(y + 1) yhat = np.log(yhat + 1) @@ -2436,12 +3270,12 @@ def plot_heatmap( ax = fig.add_subplot(111) # Build heatmap matrix. - ## Add in data. + # Add in data. data = np.array([ - self._continuous_model(idx=g, nonnumeric=False) + self._continuous_model(idx=g, non_numeric=False) for i, g in enumerate(gene_idx) ]) - ## Order columns by continuous covariate. + # Order columns by continuous covariate. idx_x_sorted = np.argsort(self._continuous_coords) data = data[:, idx_x_sorted] xcoord = self._continuous_coords[idx_x_sorted] @@ -2487,7 +3321,7 @@ def plot_heatmap( plt.close(fig) if return_axs: - return axs + return ax else: return @@ -2500,30 +3334,35 @@ def __init__( de_test: DifferentialExpressionTestWald, size_factors: np.ndarray, continuous_coords: np.ndarray, - spline_coefs: list + spline_coefs: list, + noise_model: str, ): super().__init__( de_test=de_test, model_estim=de_test.model_estim, size_factors=size_factors, continuous_coords=continuous_coords, - spline_coefs=spline_coefs + spline_coefs=spline_coefs, + noise_model=noise_model ) class DifferentialExpressionTestLRTCont(_DifferentialExpressionTestCont): de_test: DifferentialExpressionTestLRT - def __init__(self, + def __init__( + self, de_test: DifferentialExpressionTestLRT, size_factors: np.ndarray, continuous_coords: np.ndarray, - spline_coefs: list + spline_coefs: list, + noise_model: str ): super().__init__( de_test=de_test, model_estim=de_test.full_estim, size_factors=size_factors, continuous_coords=continuous_coords, - spline_coefs=spline_coefs + spline_coefs=spline_coefs, + noise_model=noise_model ) diff --git a/diffxpy/testing/tests.py b/diffxpy/testing/tests.py index 1c0beda..cba2037 100644 --- a/diffxpy/testing/tests.py +++ b/diffxpy/testing/tests.py @@ -8,6 +8,11 @@ import scipy.sparse import xarray as xr +try: + from anndata.base import Raw +except ImportError: + from anndata import Raw + from batchglm import data as data_utils from batchglm.xarray_sparse import SparseXArrayDataSet from diffxpy import pkg_constants @@ -17,9 +22,8 @@ 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") +from .utils import parse_gene_names, parse_data, parse_sample_description, parse_size_factors, parse_grouping, \ + constraint_system_from_star # Use this to suppress matrix subclass PendingDepreceationWarnings from numpy: np.warnings.filterwarnings("ignore") @@ -145,8 +149,8 @@ def _fit( else: raise ValueError('base.test(): `noise_model="%s"` not recognized.' % noise_model) - logger.info("Fitting model...") - logger.debug(" * Assembling input data...") + logging.getLogger("diffxpy").info("Fitting model...") + logging.getLogger("diffxpy").debug(" * Assembling input data...") input_data = InputData.new( data=data, design_loc=design_loc, @@ -157,7 +161,7 @@ def _fit( feature_names=gene_names, ) - logger.debug(" * Set up Estimator...") + logging.getLogger("diffxpy").debug(" * Set up Estimator...") constructor_args = {} if batch_size is not None: constructor_args["batch_size"] = batch_size @@ -176,10 +180,10 @@ def _fit( **constructor_args ) - logger.debug(" * Initializing Estimator...") + logging.getLogger("diffxpy").debug(" * Initializing Estimator...") estim.initialize() - logger.debug(" * Run estimation...") + logging.getLogger("diffxpy").debug(" * Run estimation...") # training: if callable(training_strategy): # call training_strategy if it is a function @@ -188,28 +192,28 @@ def _fit( estim.train_sequence(training_strategy=training_strategy) if close_session: - logger.debug(" * Finalize estimation...") + logging.getLogger("diffxpy").debug(" * Finalize estimation...") model = estim.finalize() else: model = estim - logger.debug(" * Model fitting done.") + logging.getLogger("diffxpy").debug(" * Model fitting done.") return model def lrt( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], - reduced_formula_loc: str, + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], full_formula_loc: str, - reduced_formula_scale: str = "~1", + reduced_formula_loc: str, full_formula_scale: str = "~1", + reduced_formula_scale: str = "~1", as_numeric: Union[List[str], Tuple[str], str] = (), init_a: Union[np.ndarray, str] = "AUTO", init_b: Union[np.ndarray, str] = "AUTO", - gene_names=None, + gene_names: Union[np.ndarray, list] = None, sample_description: pd.DataFrame = None, noise_model="nb", - size_factors: np.ndarray = None, + size_factors: Union[np.ndarray, pd.core.series.Series, np.ndarray] = None, batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "DEFAULT", quick_scale: bool = False, @@ -222,20 +226,19 @@ def lrt( Note that lrt() does not support constraints in its current form. Please use wald() for constraints. - :param data: Array-like, xr.DataArray, xr.Dataset or anndata.Anndata object containing observations. - Input data matrix (observations x features) or (cells x genes). - :param reduced_formula_loc: formula - Reduced model formula for location and scale parameter models. - If not specified, `reduced_formula` will be used instead. + :param data: Input data matrix (observations x features) or (cells x genes). :param full_formula_loc: formula Full model formula for location parameter model. If not specified, `full_formula` will be used instead. - :param reduced_formula_scale: formula - Reduced model formula for scale parameter model. + :param reduced_formula_loc: formula + Reduced model formula for location and scale parameter models. If not specified, `reduced_formula` will be used instead. :param full_formula_scale: formula Full model formula for scale parameter model. If not specified, `reduced_formula_scale` will be used instead. + :param reduced_formula_scale: formula + Reduced model formula for scale parameter model. + If not specified, `reduced_formula` will be used instead. :param as_numeric: Which columns of sample_description to treat as numeric and not as categorical. This yields columns in the design matrix @@ -247,7 +250,6 @@ def lrt( - str: * "auto": automatically choose best initialization - * "random": initialize with random values * "standard": initialize intercept with observed mean * "init_model": initialize with another model (see `ìnit_model` parameter) * "closed_form": try to initialize with closed form @@ -257,7 +259,6 @@ def lrt( - str: * "auto": automatically choose best initialization - * "random": initialize with random values * "standard": initialize with zeros * "init_model": initialize with another model (see `ìnit_model` parameter) * "closed_form": try to initialize with closed form @@ -268,7 +269,8 @@ def lrt( - 'nb': default :param size_factors: 1D array of transformed library size factors for each cell in the - same order as in data + same order as in data or string-type column identifier of size-factor containing + column in sample description. :param batch_size: the batch size to use for the estimator :param training_strategy: {str, function, list} training strategy to use. Can be: @@ -298,7 +300,7 @@ def lrt( """ # TODO test nestedness if len(kwargs) != 0: - logger.info("additional kwargs: %s", str(kwargs)) + logging.getLogger("diffxpy").info("additional kwargs: %s", str(kwargs)) if isinstance(as_numeric, str): as_numeric = [as_numeric] @@ -306,27 +308,35 @@ def lrt( 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) + size_factors = parse_size_factors( + size_factors=size_factors, + data=X, + sample_description=sample_description + ) full_design_loc = data_utils.design_matrix( sample_description=sample_description, formula=full_formula_loc, - as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] + as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values], + return_type="patsy" ) reduced_design_loc = data_utils.design_matrix( sample_description=sample_description, formula=reduced_formula_loc, - as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] + as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values], + return_type="patsy" ) full_design_scale = data_utils.design_matrix( sample_description=sample_description, formula=full_formula_scale, - as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] + as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values], + return_type="patsy" ) reduced_design_scale = data_utils.design_matrix( sample_description=sample_description, formula=reduced_formula_scale, - as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] + as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values], + return_type="patsy" ) reduced_model = _fit( @@ -377,22 +387,22 @@ def lrt( def wald( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], factor_loc_totest: Union[str, List[str]] = None, coef_to_test: Union[str, List[str]] = None, - formula_loc: str = None, - formula_scale: str = "~1", + formula_loc: Union[None, str] = None, + formula_scale: Union[None, str] = "~1", as_numeric: Union[List[str], Tuple[str], str] = (), init_a: Union[np.ndarray, str] = "AUTO", init_b: Union[np.ndarray, str] = "AUTO", - gene_names: Union[str, np.ndarray] = None, - sample_description: pd.DataFrame = None, + gene_names: Union[np.ndarray, list] = None, + sample_description: Union[None, pd.DataFrame] = None, dmat_loc: Union[patsy.design_info.DesignMatrix, xr.Dataset] = None, dmat_scale: Union[patsy.design_info.DesignMatrix, xr.Dataset] = None, - constraints_loc: np.ndarray = None, - constraints_scale: np.ndarray = None, + constraints_loc: Union[None, List[str], Tuple[str, str], dict, np.ndarray] = None, + constraints_scale: Union[None, List[str], Tuple[str, str], dict, np.ndarray] = None, noise_model: str = "nb", - size_factors: np.ndarray = None, + size_factors: Union[np.ndarray, pd.core.series.Series, str] = None, batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO", quick_scale: bool = False, @@ -402,8 +412,7 @@ def wald( """ Perform Wald test for differential expression for each gene. - :param data: Array-like, xr.DataArray, xr.Dataset or anndata.Anndata object containing observations. - Input data matrix (observations x features) or (cells x genes). + :param data: Input data matrix (observations x features) or (cells x genes). :param factor_loc_totest: str, list of strings List of factors of formula to test with Wald test. E.g. "condition" or ["batch", "condition"] if formula_loc would be "~ 1 + batch + condition" @@ -412,8 +421,6 @@ def wald( this parameter allows to specify the group which should be tested. Alternatively, if factor_loc_totest is not given, this list sets the exact coefficients which are to be tested. - :param formula: formula - model formula for location and scale parameter models. :param formula_loc: formula model formula for location and scale parameter models. If not specified, `formula` will be used instead. @@ -423,7 +430,7 @@ def wald( :param as_numeric: Which columns of sample_description to treat as numeric and not as categorical. This yields columns in the design matrix - which do not correpond to one-hot encoded discrete factors. + which do not correspond to one-hot encoded discrete factors. This makes sense for number of genes, time, pseudotime or space for example. :param init_a: (Optional) Low-level initial values for a. @@ -450,32 +457,79 @@ def wald( :param dmat_scale: Pre-built scale model design matrix. This over-rides formula_scale and sample description information given in data or sample_description. - :param constraints_loc: : Constraints for location model. - Array with constraints in rows and model parameters in columns. - Each constraint contains non-zero entries for the a of parameters that - has to sum to zero. This constraint is enforced by binding one parameter - to the negative sum of the other parameters, effectively representing that - parameter as a function of the other parameters. This dependent - parameter is indicated by a -1 in this array, the independent parameters - of that constraint (which may be dependent at an earlier constraint) - are indicated by a 1. It is highly recommended to only use this option - together with prebuilt design matrix for the location model, dmat_loc. - :param constraints_scale: : Constraints for scale model. - Array with constraints in rows and model parameters in columns. - Each constraint contains non-zero entries for the a of parameters that - has to sum to zero. This constraint is enforced by binding one parameter - to the negative sum of the other parameters, effectively representing that - parameter as a function of the other parameters. This dependent - parameter is indicated by a -1 in this array, the independent parameters - of that constraint (which may be dependent at an earlier constraint) - are indicated by a 1. It is highly recommended to only use this option - together with prebuilt design matrix for the scale model, dmat_scale. + :param constraints_loc: Constraints for location model. Can be one of the following: + + - np.ndarray: + Array with constraints in rows and model parameters in columns. + Each constraint contains non-zero entries for the a of parameters that + has to sum to zero. This constraint is enforced by binding one parameter + to the negative sum of the other parameters, effectively representing that + parameter as a function of the other parameters. This dependent + parameter is indicated by a -1 in this array, the independent parameters + of that constraint (which may be dependent at an earlier constraint) + are indicated by a 1. You should only use this option + together with prebuilt design matrix for the location model, dmat_loc, + for example via de.utils.setup_constrained(). + - dict: + Every element of the dictionary corresponds to one set of equality constraints. + Each set has to be be an entry of the form {..., x: y, ...} + where x is the factor to be constrained and y is a factor by which levels of x are grouped + and then constrained. Set y="1" to constrain all levels of x to sum to one, + a single equality constraint. + + E.g.: {"batch": "condition"} Batch levels within each condition are constrained to sum to + zero. This is applicable if repeats of a an experiment within each condition + are independent so that the set-up ~1+condition+batch is perfectly confounded. + + Can only group by non-constrained effects right now, use constraint_matrix_from_string + for other cases. + - list of strings or tuple of strings: + String encoded equality constraints. + + E.g. ["batch1 + batch2 + batch3 = 0"] + - None: + No constraints are used, this is equivalent to using an identity matrix as a + constraint matrix. + :param constraints_scale: Constraints for scale model. Can be one of the following: + + - np.ndarray: + Array with constraints in rows and model parameters in columns. + Each constraint contains non-zero entries for the a of parameters that + has to sum to zero. This constraint is enforced by binding one parameter + to the negative sum of the other parameters, effectively representing that + parameter as a function of the other parameters. This dependent + parameter is indicated by a -1 in this array, the independent parameters + of that constraint (which may be dependent at an earlier constraint) + are indicated by a 1. You should only use this option + together with prebuilt design matrix for the scale model, dmat_scale, + for example via de.utils.setup_constrained(). + - dict: + Every element of the dictionary corresponds to one set of equality constraints. + Each set has to be be an entry of the form {..., x: y, ...} + where x is the factor to be constrained and y is a factor by which levels of x are grouped + and then constrained. Set y="1" to constrain all levels of x to sum to one, + a single equality constraint. + + E.g.: {"batch": "condition"} Batch levels within each condition are constrained to sum to + zero. This is applicable if repeats of a an experiment within each condition + are independent so that the set-up ~1+condition+batch is perfectly confounded. + + Can only group by non-constrained effects right now, use constraint_matrix_from_string + for other cases. + - list of strings or tuple of strings: + String encoded equality constraints. + + E.g. ["batch1 + batch2 + batch3 = 0"] + - None: + No constraints are used, this is equivalent to using an identity matrix as a + constraint matrix. :param size_factors: 1D array of transformed library size factors for each cell in the - same order as in data + same order as in data or string-type column identifier of size-factor containing + column in sample description. :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default - :param batch_size: the batch size to use for the estimator + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -483,17 +537,6 @@ def wald( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. :param quick_scale: Depending on the optimizer, `scale` will be fitted faster and maybe less accurate. Useful in scenarios where fitting the exact `scale` is not absolutely necessary. @@ -503,12 +546,16 @@ def wald( :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ if len(kwargs) != 0: - logger.debug("additional kwargs: %s", str(kwargs)) - - if dmat_loc is None and formula_loc is None: - raise ValueError("Supply either dmat_loc or formula_loc or formula.") - if dmat_scale is None and formula_scale is None: - raise ValueError("Supply either dmat_loc or formula_loc or formula.") + logging.getLogger("diffxpy").debug("additional kwargs: %s", str(kwargs)) + + if (dmat_loc is None and formula_loc is None) or \ + (dmat_loc is not None and formula_loc is not None): + raise ValueError("Supply either dmat_loc or formula_loc.") + if (dmat_scale is None and formula_scale is None) or \ + (dmat_scale is not None and formula_scale != "~1"): + raise ValueError("Supply either dmat_scale or formula_scale.") + if dmat_loc is not None and factor_loc_totest is not None: + raise ValueError("Supply coef_to_test and not factor_loc_totest if dmat_loc is supplied.") # Check that factor_loc_totest and coef_to_test are lists and not single strings: if isinstance(factor_loc_totest, str): factor_loc_totest = [factor_loc_totest] @@ -522,46 +569,32 @@ def wald( 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) - - if dmat_loc is None: - design_loc = data_utils.design_matrix( - sample_description=sample_description, - formula=formula_loc, - as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] - ) - # Check that closed-form is not used if numeric predictors are used and model is not "norm". - if isinstance(init_a, str): - if np.any([True if x in as_numeric else False for x in sample_description.columns.values]): - if noise_model.lower() not in ["normal", "norm"]: - if init_a == "closed_form": - init_a = "standard" - logger.warning("Setting init_a to standard as numeric predictors were supplied.") - logger.warning("Closed-form initialisation is not possible" + - " for noise model %s with numeric predictors." % noise_model) - elif init_a == "AUTO": - init_a = "standard" - else: - design_loc = dmat_loc + size_factors = parse_size_factors( + size_factors=size_factors, + data=X, + sample_description=sample_description + ) - if dmat_scale is None: - design_scale = data_utils.design_matrix( - sample_description=sample_description, - formula=formula_scale, - as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] - ) - # Check that closed-form is not used if numeric predictors are used and model is not "norm". - if isinstance(init_b, str): - if np.any([True if x in as_numeric else False for x in sample_description.columns.values]): - if init_b == "closed_form": - init_b = "standard" - logger.warning("Setting init_b to standard as numeric predictors were supplied.") - logger.warning("Closed-form initialisation is not possible" + - " for noise model %s with numeric predictors." % noise_model) - elif init_b == "AUTO": - init_b = "standard" - else: - design_scale = dmat_scale + logging.getLogger("diffxpy").debug("building location model") + design_loc, constraints_loc = constraint_system_from_star( + dmat=dmat_loc, + sample_description=sample_description, + formula=formula_loc, + as_numeric=as_numeric, + constraints=constraints_loc, + dims=["design_loc_params", "loc_params"], + return_type="patsy" + ) + logging.getLogger("diffxpy").debug("building scale model") + design_scale, constraints_scale = constraint_system_from_star( + dmat=dmat_scale, + sample_description=sample_description, + formula=formula_scale, + as_numeric=as_numeric, + constraints=constraints_scale, + dims=["design_scale_params", "scale_params"], + return_type="patsy" + ) # Define indices of coefficients to test: constraints_loc_temp = constraints_loc if constraints_loc is not None else np.eye(design_loc.shape[-1]) @@ -583,18 +616,16 @@ def wald( elif coef_to_test is not None: # Directly select coefficients to test from design matrix (xarray): # Check that coefficients to test are not dependent parameters if constraints are given: - # TODO: design_loc is sometimes xarray and sometimes patsy when it arrives here, - # should it not always be xarray? - if isinstance(design_loc, patsy.design_info.DesignMatrix): - col_indices = np.asarray([ - design_loc.design_info.column_names.index(x) - for x in coef_to_test - ]) - else: - col_indices = np.asarray([ - list(np.asarray(design_loc.coords['design_params'])).index(x) - for x in coef_to_test - ]) + coef_loc_names = data_utils.view_coef_names(design_loc).tolist() + if not np.all([x in coef_loc_names for x in coef_to_test]): + raise ValueError( + "the requested test coefficients %s were found in model coefficients %s" % + (", ".join([x for x in coef_to_test if x not in coef_loc_names]), + ", ".join(coef_loc_names)) + ) + col_indices = np.asarray([ + coef_loc_names.index(x) for x in coef_to_test + ]) else: raise ValueError("either set factor_loc_totest or coef_to_test") # Check that all tested coefficients are independent: @@ -624,18 +655,21 @@ def wald( de_test = DifferentialExpressionTestWald( model_estim=model, - col_indices=col_indices + col_indices=col_indices, + noise_model=noise_model, + sample_description=sample_description ) return de_test def t_test( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], grouping, - gene_names=None, - sample_description=None, - is_logged=False, + gene_names: Union[np.ndarray, list] = None, + sample_description: pd.DataFrame = None, + is_logged: bool = False, + is_sig_zerovar: bool = True, dtype="float64" ): """ @@ -653,6 +687,9 @@ def t_test( :param is_logged: 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. + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. """ gene_names = parse_gene_names(data, gene_names) X = parse_data(data, gene_names) @@ -662,20 +699,23 @@ def t_test( de_test = DifferentialExpressionTestTT( data=X.astype(dtype), + sample_description=sample_description, grouping=grouping, gene_names=gene_names, - is_logged=is_logged + is_logged=is_logged, + is_sig_zerovar=is_sig_zerovar ) return de_test def rank_test( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], - grouping, - gene_names=None, - sample_description=None, - is_logged=False, + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + grouping: Union[str, np.ndarray, list], + gene_names: Union[np.ndarray, list] = None, + sample_description: pd.DataFrame = None, + is_logged: bool = False, + is_sig_zerovar: bool = True, dtype="float64" ): """ @@ -693,6 +733,9 @@ def rank_test( :param is_logged: 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. + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. """ gene_names = parse_gene_names(data, gene_names) X = parse_data(data, gene_names) @@ -702,25 +745,28 @@ def rank_test( de_test = DifferentialExpressionTestRank( data=X.astype(dtype), + sample_description=sample_description, grouping=grouping, gene_names=gene_names, - is_logged=is_logged + is_logged=is_logged, + is_sig_zerovar=is_sig_zerovar ) return de_test def two_sample( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], grouping: Union[str, np.ndarray, list], as_numeric: Union[List[str], Tuple[str], str] = (), - test=None, - gene_names=None, - sample_description=None, + test: str = "t-test", + gene_names: Union[np.ndarray, list] = None, + sample_description: pd.DataFrame = None, noise_model: str = None, size_factors: np.ndarray = None, batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO", + is_sig_zerovar: bool = True, quick_scale: bool = None, dtype="float64", **kwargs @@ -735,7 +781,7 @@ def two_sample( The exact unit_test are as follows (assuming the group labels are saved in a column named "group"): - - lrt(log-likelihood ratio test): + - "lrt" - (log-likelihood ratio test): Requires the fitting of 2 generalized linear models (full and reduced). The models are automatically assembled as follows, use the de.test.lrt() function if you would like to perform a different test. @@ -744,15 +790,15 @@ def two_sample( * full model scale parameter: ~ 1 + group * reduced model location parameter: ~ 1 * reduced model scale parameter: ~ 1 + group - - Wald test: + - "wald" - Wald test: Requires the fitting of 1 generalized linear models. model location parameter: ~ 1 + group model scale parameter: ~ 1 + group Test the group coefficient of the location parameter model against 0. - - t-test: + - "t-test" - Welch's t-test: Doesn't require fitting of generalized linear models. Welch's t-test between both observation groups. - - wilcoxon: + - "rank" - Wilcoxon rank sum (Mann-Whitney U) test: Doesn't require fitting of generalized linear models. Wilcoxon rank sum (Mann-Whitney U) test between both observation groups. @@ -773,15 +819,15 @@ def two_sample( - 'wald': default - 'lrt' - 't-test' - - 'wilcoxon' + - 'rank' :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 + :param size_factors: 1D array of transformed library size factors for each cell in the + same order as in data :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default - :param size_factors: 1D array of transformed library size factors for each cell in the - same order as in data - :param batch_size: the batch size to use for the estimator + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -789,17 +835,9 @@ def two_sample( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. :param quick_scale: Depending on the optimizer, `scale` will be fitted faster and maybe less accurate. Useful in scenarios where fitting the exact `scale` is not absolutely necessary. @@ -808,8 +846,8 @@ def two_sample( Should be "float32" for single precision or "float64" for double precision. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ - if test in ['t-test', 'wilcoxon'] and noise_model is not None: - raise ValueError('base.two_sample(): Do not specify `noise_model` if using test t-test or wilcoxon: ' + + if test in ['t-test', 'rank'] and noise_model is not None: + raise ValueError('base.two_sample(): Do not specify `noise_model` if using test t-test or rank_test: ' + 'The t-test is based on a gaussian noise model and wilcoxon is model free.') gene_names = parse_gene_names(data, gene_names) @@ -823,10 +861,6 @@ def two_sample( if groups.size < 2: raise ValueError("Less than two groups detected:\n\t%s", groups) - # Set default test: - if test is None: - test = 'wald' - if test.lower() == 'wald': if noise_model is None: raise ValueError("Please specify noise_model") @@ -853,9 +887,9 @@ def two_sample( if noise_model is None: raise ValueError("Please specify noise_model") full_formula_loc = '~ 1 + grouping' - full_formula_scale = '~ 1 + grouping' + full_formula_scale = '~ 1' reduced_formula_loc = '~ 1' - reduced_formula_scale = '~ 1 + grouping' + reduced_formula_scale = '~ 1' de_test = lrt( data=X, full_formula_loc=full_formula_loc, @@ -878,36 +912,39 @@ def two_sample( data=X, gene_names=gene_names, grouping=grouping, + is_sig_zerovar=is_sig_zerovar, dtype=dtype ) - elif test.lower() == 'wilcoxon': + elif test.lower() == 'rank': de_test = rank_test( data=X, gene_names=gene_names, grouping=grouping, + is_sig_zerovar=is_sig_zerovar, dtype=dtype ) else: - raise ValueError('base.two_sample(): Parameter `test="%s"` not recognized.' % test) + raise ValueError('two_sample(): Parameter `test="%s"` not recognized.' % test) return de_test def pairwise( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], grouping: Union[str, np.ndarray, list], - as_numeric: Union[List[str], Tuple[str], str] = [], + as_numeric: Union[List[str], Tuple[str], str] = (), test: str = 'z-test', lazy: bool = False, - gene_names: str = None, + gene_names: Union[np.ndarray, list] = None, sample_description: pd.DataFrame = None, noise_model: str = None, - pval_correction: str = "global", size_factors: np.ndarray = None, batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO", + is_sig_zerovar: bool = True, quick_scale: bool = None, dtype="float64", + pval_correction: str = "global", keep_full_test_objs: bool = False, **kwargs ): @@ -926,22 +963,22 @@ def pairwise( on the subset of the data that only contains observations of a given pair of groups: - - lrt(log-likelihood ratio test): + - "lrt" -log-likelihood ratio test: Requires the fitting of 2 generalized linear models (full and reduced). * full model location parameter: ~ 1 + group * full model scale parameter: ~ 1 + group * reduced model location parameter: ~ 1 * reduced model scale parameter: ~ 1 + group - - Wald test: + - "wald" - Wald test: Requires the fitting of 1 generalized linear models. model location parameter: ~ 1 + group model scale parameter: ~ 1 + group Test the group coefficient of the location parameter model against 0. - - t-test: + - "t-test" - Welch's t-test: Doesn't require fitting of generalized linear models. Welch's t-test between both observation groups. - - wilcoxon: + - "rank" - Wilcoxon rank sum (Mann-Whitney U) test: Doesn't require fitting of generalized linear models. Wilcoxon rank sum (Mann-Whitney U) test between both observation groups. @@ -963,7 +1000,7 @@ def pairwise( - 'wald' - 'lrt' - 't-test' - - 'wilcoxon' + - 'rank' :param lazy: bool, whether to enable lazy results evaluation. This is only possible if test=="ztest" and yields an output object which computes p-values etc. only upon request of certain pairs. This makes sense if the entire @@ -972,17 +1009,12 @@ def pairwise( a certain subset of the pairwise comparisons is desired anyway. :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 + :param size_factors: 1D array of transformed library size factors for each cell in the + same order as in data :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default - :param pval_correction: Choose between global and test-wise correction. - Can be: - - - "global": correct all p-values in one operation - - "by_test": correct the p-values of each test individually - :param size_factors: 1D array of transformed library size factors for each cell in the - same order as in data - :param batch_size: the batch size to use for the estimator + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -990,28 +1022,25 @@ def pairwise( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. :param quick_scale: Depending on the optimizer, `scale` will be fitted faster and maybe less accurate. Useful in scenarios where fitting the exact `scale` is not absolutely necessary. :param dtype: Allows specifying the precision which should be used to fit data. Should be "float32" for single precision or "float64" for double precision. - :param keep_full_test_objs: [Debugging] keep the individual test objects; currently valid for test != "z-test" + :param pval_correction: Choose between global and test-wise correction. + Can be: + + - "global": correct all p-values in one operation + - "by_test": correct the p-values of each test individually + :param keep_full_test_objs: [Debugging] keep the individual test objects; currently valid for test != "z-test". + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ if len(kwargs) != 0: - logger.info("additional kwargs: %s", str(kwargs)) + logging.getLogger("diffxpy").info("additional kwargs: %s", str(kwargs)) if lazy and not (test.lower() == 'z-test' or test.lower() == 'z_test' or test.lower() == 'ztest'): raise ValueError("lazy evaluation of pairwise tests only possible if test is z-test") @@ -1087,6 +1116,7 @@ def pairwise( batch_size=batch_size, training_strategy=training_strategy, quick_scale=quick_scale, + is_sig_zerovar=is_sig_zerovar, dtype=dtype, **kwargs ) @@ -1112,19 +1142,20 @@ def pairwise( def versus_rest( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], grouping: Union[str, np.ndarray, list], as_numeric: Union[List[str], Tuple[str], str] = (), test: str = 'wald', - gene_names: str = None, + gene_names: Union[np.ndarray, list] = None, sample_description: pd.DataFrame = None, noise_model: str = None, - pval_correction: str = "global", size_factors: np.ndarray = None, batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO", + is_sig_zerovar: bool = True, quick_scale: bool = None, dtype="float64", + pval_correction: str = "global", keep_full_test_objs: bool = False, **kwargs ): @@ -1144,22 +1175,22 @@ def versus_rest( is one group and the remaining groups are allocated to the second reference group): - - lrt(log-likelihood ratio test): + - "lrt" - log-likelihood ratio test): Requires the fitting of 2 generalized linear models (full and reduced). * full model location parameter: ~ 1 + group * full model scale parameter: ~ 1 + group * reduced model location parameter: ~ 1 * reduced model scale parameter: ~ 1 + group - - Wald test: + - "wald" - Wald test: Requires the fitting of 1 generalized linear models. model location parameter: ~ 1 + group model scale parameter: ~ 1 + group Test the group coefficient of the location parameter model against 0. - - t-test: + - "t-test" - Welch's t-test: Doesn't require fitting of generalized linear models. Welch's t-test between both observation groups. - - wilcoxon: + - "rank" - Wilcoxon rank sum (Mann-Whitney U) test: Doesn't require fitting of generalized linear models. Wilcoxon rank sum (Mann-Whitney U) test between both observation groups. @@ -1175,17 +1206,14 @@ def versus_rest( which do not correpond to one-hot encoded discrete factors. This makes sense for number of genes, time, pseudotime or space for example. - :param test: str, statistical test to use. Possible options: + :param test: str, statistical test to use. Possible options (see function description): - 'wald' - 'lrt' - 't-test' - - 'wilcoxon' + - 'rank' :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 - :param noise_model: str, noise model to use in model-based unit_test. Possible options: - - - 'nb': default :param pval_correction: Choose between global and test-wise correction. Can be: @@ -1193,7 +1221,10 @@ def versus_rest( - "by_test": correct the p-values of each test individually :param size_factors: 1D array of transformed library size factors for each cell in the same order as in data - :param batch_size: the batch size to use for the estimator + :param noise_model: str, noise model to use in model-based unit_test. Possible options: + + - 'nb': default + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -1201,28 +1232,24 @@ def versus_rest( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. :param quick_scale: Depending on the optimizer, `scale` will be fitted faster and maybe less accurate. - Useful in scenarios where fitting the exact `scale` is not + Useful in scenarios where fitting the exact `scale` is not absolutely necessary. :param dtype: Allows specifying the precision which should be used to fit data. Should be "float32" for single precision or "float64" for double precision. - :param keep_full_test_objs: [Debugging] keep the individual test objects; currently valid for test != "z-test" + :param pval_correction: Choose between global and test-wise correction. + Can be: + + - "global": correct all p-values in one operation + - "by_test": correct the p-values of each test individually + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ if len(kwargs) != 0: - logger.info("additional kwargs: %s", str(kwargs)) + logging.getLogger("diffxpy").info("additional kwargs: %s", str(kwargs)) # Do not store all models but only p-value and q-value matrix: # genes x groups @@ -1255,6 +1282,7 @@ def versus_rest( training_strategy=training_strategy, quick_scale=quick_scale, size_factors=size_factors, + is_sig_zerovar=is_sig_zerovar, dtype=dtype, **kwargs ) @@ -1277,10 +1305,11 @@ def versus_rest( def partition( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], - partition: Union[str, np.ndarray, list], - gene_names: str = None, - sample_description: pd.DataFrame = None): + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + parts: Union[str, np.ndarray, list], + gene_names: Union[np.ndarray, list] = None, + sample_description: pd.DataFrame = None +): r""" Perform differential expression test for each group. This class handles the partitioning of the data set, the differential test callls and @@ -1293,18 +1322,21 @@ def partition( :param data: Array-like, xr.DataArray, xr.Dataset or anndata.Anndata object containing observations. Input data matrix (observations x features) or (cells x genes). + :param parts: str, array + - column in data.obs/sample_description which contains the split of observations into the two groups. + - array of length `num_observations` containing group labels :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 """ return (_Partition( data=data, - partition=partition, + parts=parts, gene_names=gene_names, sample_description=sample_description)) -class _Partition(): +class _Partition: """ Perform differential expression test for each group. This class handles the partitioning of the data set, the differential test callls and @@ -1317,13 +1349,14 @@ class _Partition(): def __init__( self, data: Union[anndata.AnnData, xr.DataArray, xr.Dataset, np.ndarray], - partition: Union[str, np.ndarray, list], - gene_names: str = None, - sample_description: pd.DataFrame = None): + parts: Union[str, np.ndarray, list], + gene_names: Union[np.ndarray, list] = None, + sample_description: pd.DataFrame = None + ): """ :param data: Array-like, xr.DataArray, xr.Dataset or anndata.Anndata object containing observations. Input data matrix (observations x features) or (cells x genes). - :param partition: str, array + :param parts: str, array - column in data.obs/sample_description which contains the split of observations into the two groups. - array of length `num_observations` containing group labels @@ -1333,7 +1366,7 @@ def __init__( 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.partition = parse_grouping(data, sample_description, parts) self.partitions = np.unique(self.partition) self.partition_idx = [np.where(self.partition == x)[0] for x in self.partitions] @@ -1342,14 +1375,17 @@ def two_sample( grouping: Union[str], as_numeric: Union[List[str], Tuple[str], str] = (), test=None, - noise_model: str = None, size_factors: np.ndarray = None, + noise_model: str = None, batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO", + is_sig_zerovar: bool = True, **kwargs ) -> _DifferentialExpressionTestMulti: """ - See annotation of de.test.two_sample() + Performs a two-sample test within each partition of a data set. + + See also annotation of de.test.two_sample() :param grouping: str @@ -1365,11 +1401,13 @@ def two_sample( - 'wald': default - 'lrt' - 't-test' - - 'wilcoxon' + - 'rank' + :param size_factors: 1D array of transformed library size factors for each cell in the + same order as in data :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default - :param batch_size: the batch size to use for the estimator + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -1377,17 +1415,9 @@ def two_sample( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ DETestsSingle = [] @@ -1403,6 +1433,7 @@ def two_sample( size_factors=size_factors[idx] if size_factors is not None else None, batch_size=batch_size, training_strategy=training_strategy, + is_sig_zerovar=is_sig_zerovar, **kwargs )) return DifferentialExpressionTestByPartition( @@ -1414,22 +1445,35 @@ def two_sample( def t_test( self, grouping: Union[str], + is_logged: bool, + is_sig_zerovar: bool = True, dtype="float64" ): """ - See annotation of de.test.t_test() + Performs a Welch's t-test within each partition of a data set. + + See also annotation of de.test.t_test() :param grouping: str - column in data.obs/sample_description which contains the split of observations into the two groups. + :param is_logged: + 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. + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. + :param dtype: """ DETestsSingle = [] for i, idx in enumerate(self.partition_idx): DETestsSingle.append(t_test( data=self.X[idx, :], grouping=grouping, + is_logged=is_logged, gene_names=self.gene_names, sample_description=self.sample_description.iloc[idx, :], + is_sig_zerovar=is_sig_zerovar, dtype=dtype )) return DifferentialExpressionTestByPartition( @@ -1438,18 +1482,25 @@ def t_test( ave=np.mean(self.X, axis=0), correction_type="by_test") - def wilcoxon( + def rank_test( self, grouping: Union[str], + is_sig_zerovar: bool = True, dtype="float64" ): """ - See annotation of de.test.wilcoxon() + Performs a Wilcoxon rank sum test within each partition of a data set. + + See also annotation of de.test.rank_test() :param grouping: str, array - column in data.obs/sample_description which contains the split of observations into the two groups. - array of length `num_observations` containing group labels + :param is_sig_zerovar: + Whether to assign p-value of 0 to a gene which has zero variance in both groups but not the same mean. If False, + the p-value is set to np.nan. + :param dtype: """ DETestsSingle = [] for i, idx in enumerate(self.partition_idx): @@ -1458,6 +1509,7 @@ def wilcoxon( grouping=grouping, gene_names=self.gene_names, sample_description=self.sample_description.iloc[idx, :], + is_sig_zerovar=is_sig_zerovar, dtype=dtype )) return DifferentialExpressionTestByPartition( @@ -1468,44 +1520,68 @@ def wilcoxon( def lrt( self, - reduced_formula_loc: str = None, full_formula_loc: str = None, - reduced_formula_scale: str = None, + reduced_formula_loc: str = None, full_formula_scale: str = "~1", + reduced_formula_scale: str = None, as_numeric: Union[List[str], Tuple[str], str] = (), - noise_model="nb", + init_a: Union[str] = "AUTO", + init_b: Union[str] = "AUTO", size_factors: np.ndarray = None, + noise_model="nb", batch_size: int = None, training_strategy: Union[str, List[Dict[str, object]], Callable] = "AUTO", **kwargs ): """ - See annotation of de.test.lrt() + Performs a likelihood-ratio test within each partition of a data set. + + See also annotation of de.test.lrt() - :param reduced_formula_loc: formula - Reduced model formula for location and scale parameter models. - If not specified, `reduced_formula` will be used instead. :param full_formula_loc: formula Full model formula for location parameter model. If not specified, `full_formula` will be used instead. - :param reduced_formula_scale: formula - Reduced model formula for scale parameter model. + :param reduced_formula_loc: formula + Reduced model formula for location and scale parameter models. If not specified, `reduced_formula` will be used instead. :param full_formula_scale: formula Full model formula for scale parameter model. If not specified, `reduced_formula_scale` will be used instead. + :param reduced_formula_scale: formula + Reduced model formula for scale parameter model. + If not specified, `reduced_formula` will be used instead. :param as_numeric: Which columns of sample_description to treat as numeric and not as categorical. This yields columns in the design matrix which do not correpond to one-hot encoded discrete factors. This makes sense for number of genes, time, pseudotime or space for example. + :param init_a: (Optional) Low-level initial values for a. + Can be: + + - str: + * "auto": automatically choose best initialization + * "standard": initialize intercept with observed mean + * "init_model": initialize with another model (see `ìnit_model` parameter) + * "closed_form": try to initialize with closed form + + Note that unlike in the lrt without partitions, this does not support np.ndarrays. + :param init_b: (Optional) Low-level initial values for b + Can be: + + - str: + * "auto": automatically choose best initialization + * "standard": initialize with zeros + * "init_model": initialize with another model (see `ìnit_model` parameter) + * "closed_form": try to initialize with closed form + + Note that unlike in the lrt without partitions, this does not support np.ndarrays. + :param size_factors: 1D array of transformed library size factors for each cell in the + same order as in data :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default - :param size_factors: 1D array of transformed library size factors for each cell in the - same order as in data - :param batch_size: the batch size to use for the estimator + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -1513,17 +1589,6 @@ def lrt( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ DETestsSingle = [] @@ -1534,6 +1599,9 @@ def lrt( full_formula_loc=full_formula_loc, reduced_formula_scale=reduced_formula_scale, full_formula_scale=full_formula_scale, + as_numeric=as_numeric, + init_a=init_a, + init_b=init_b, gene_names=self.gene_names, sample_description=self.sample_description.iloc[idx, :], noise_model=noise_model, @@ -1555,6 +1623,8 @@ def wald( formula_loc: str = None, formula_scale: str = "~1", as_numeric: Union[List[str], Tuple[str], str] = (), + constraints_loc: np.ndarray = None, + constraints_scale: np.ndarray = None, noise_model: str = "nb", size_factors: np.ndarray = None, batch_size: int = None, @@ -1562,32 +1632,56 @@ def wald( **kwargs ): """ - This function performs a wald test within each partition of a data set. - See annotation of de.test.wald() - + Performs a wald test within each partition of a data set. + + See also annotation of de.test.wald() + + :param factor_loc_totest: str, list of strings + List of factors of formula to test with Wald test. + E.g. "condition" or ["batch", "condition"] if formula_loc would be "~ 1 + batch + condition" + :param coef_to_test: + If there are more than two groups specified by `factor_loc_totest`, + this parameter allows to specify the group which should be tested. + Alternatively, if factor_loc_totest is not given, this list sets + the exact coefficients which are to be tested. :param formula_loc: formula model formula for location and scale parameter models. If not specified, `formula` will be used instead. :param formula_scale: formula model formula for scale parameter model. If not specified, `formula` will be used instead. - :param factor_loc_totest: str - Factor of formula to test with Wald test. - E.g. "condition" if formula_loc would be "~ 1 + batch + condition" - :param coef_to_test: If there are more than two groups specified by `factor_loc_totest`, - this parameter allows to specify the group which should be tested :param as_numeric: Which columns of sample_description to treat as numeric and not as categorical. This yields columns in the design matrix which do not correpond to one-hot encoded discrete factors. This makes sense for number of genes, time, pseudotime or space for example. + :param constraints_loc: : Constraints for location model. + Array with constraints in rows and model parameters in columns. + Each constraint contains non-zero entries for the a of parameters that + has to sum to zero. This constraint is enforced by binding one parameter + to the negative sum of the other parameters, effectively representing that + parameter as a function of the other parameters. This dependent + parameter is indicated by a -1 in this array, the independent parameters + of that constraint (which may be dependent at an earlier constraint) + are indicated by a 1. It is highly recommended to only use this option + together with prebuilt design matrix for the location model, dmat_loc. + :param constraints_scale: : Constraints for scale model. + Array with constraints in rows and model parameters in columns. + Each constraint contains non-zero entries for the a of parameters that + has to sum to zero. This constraint is enforced by binding one parameter + to the negative sum of the other parameters, effectively representing that + parameter as a function of the other parameters. This dependent + parameter is indicated by a -1 in this array, the independent parameters + of that constraint (which may be dependent at an earlier constraint) + are indicated by a 1. It is highly recommended to only use this option + together with prebuilt design matrix for the scale model, dmat_scale. + :param size_factors: 1D array of transformed library size factors for each cell in the + same order as in data :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default - :param size_factors: 1D array of transformed library size factors for each cell in the - same order as in data - :param batch_size: the batch size to use for the estimator + :param batch_size: The batch size to use for the estimator. :param training_strategy: {str, function, list} training strategy to use. Can be: - str: will use Estimator.TrainingStrategy[training_strategy] to train @@ -1595,17 +1689,6 @@ def wald( `training_strategy(estimator)`. - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of method arguments. - - Example: - - .. code-block:: python - - [ - {"learning_rate": 0.5, }, - {"learning_rate": 0.05, }, - ] - - This will run training first with learning rate = 0.5 and then with learning rate = 0.05. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ DETestsSingle = [] @@ -1617,6 +1700,8 @@ def wald( formula_loc=formula_loc, formula_scale=formula_scale, as_numeric=as_numeric, + constraints_loc=constraints_loc, + constraints_scale=constraints_scale, gene_names=self.gene_names, sample_description=self.sample_description.iloc[idx, :], noise_model=noise_model, @@ -1633,7 +1718,7 @@ def wald( def continuous_1d( - data: Union[anndata.AnnData, anndata.base.Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], continuous: str, df: int = 5, factor_loc_totest: Union[str, List[str]] = None, @@ -1643,8 +1728,10 @@ def continuous_1d( test: str = 'wald', init_a: Union[np.ndarray, str] = "standard", init_b: Union[np.ndarray, str] = "standard", - gene_names=None, + gene_names: Union[np.ndarray, list] = None, sample_description=None, + constraints_loc: Union[dict, None] = None, + constraints_scale: Union[dict, None] = None, noise_model: str = 'nb', size_factors: np.ndarray = None, batch_size: int = None, @@ -1658,7 +1745,7 @@ def continuous_1d( This function wraps the selected statistical test for scenarios with continuous covariates and performs the necessary - spline basis transformation of the continuous covariate so that the + spline basis transformation of the continuous co-variate so that the problem can be framed as a GLM. Note that direct supply of dmats is not enabled as this function wraps @@ -1667,6 +1754,11 @@ def continuous_1d( perform these spline basis transforms outside of diffxpy and feed the dmat directly to one of the test routines wald() or lrt(). + The constraint interface only-supports dictionary-formatted constraints and + string-formatted constraints but not array-formatted constraint matrices as + design matrices are built within this function and the shape of constraint + matrices depends on the output of this function. + :param data: Array-like, xr.DataArray, xr.Dataset or anndata.Anndata object containing observations. Input data matrix (observations x features) or (cells x genes). :param continuous: str @@ -1679,11 +1771,6 @@ def continuous_1d( :param factor_loc_totest: List of factors of formula to test with Wald test. E.g. "condition" or ["batch", "condition"] if formula_loc would be "~ 1 + batch + condition" - :param formula: formula - Model formula for location and scale parameter models. - Refer to continuous covariate by the name givne in the parameter continuous, - this will be propagated across all coefficients which represent this covariate - in the spline basis space. :param formula_loc: formula Model formula for location and scale parameter models. If not specified, `formula` will be used instead. @@ -1697,11 +1784,10 @@ def continuous_1d( this will be propagated across all coefficients which represent this covariate in the spline basis space. :param as_numeric: - Which columns of sample_description to treat as numeric and - not as categorical. This yields columns in the design matrix - which do not correpond to one-hot encoded discrete factors. - This makes sense for number of genes, time, pseudotime or space - for example. + Which columns of sample_description to treat as numeric and not as categorical. + This yields columns in the design matrix which do not correpond to one-hot encoded discrete factors. + This makes sense for library depth for example. Do not use this for the covariate that you + want to extrpolate with using a spline-basis! :param test: str, statistical test to use. Possible options: - 'wald': default @@ -1720,8 +1806,59 @@ def continuous_1d( * "auto": automatically choose best initialization * "standard": initialize with zeros - np.ndarray: direct initialization of 'b' - :param gene_names: optional list/array of gene names which will be used if `data` does not implicitly store these + :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 + :param constraints_loc: Constraints for location model. Can be one of the following: + + - dict: + Every element of the dictionary corresponds to one set of equality constraints. + Each set has to be be an entry of the form {..., x: y, ...} + where x is the factor to be constrained and y is a factor by which levels of x are grouped + and then constrained. Set y="1" to constrain all levels of x to sum to one, + a single equality constraint. + + E.g.: {"batch": "condition"} Batch levels within each condition are constrained to sum to + zero. This is applicable if repeats of a an experiment within each condition + are independent so that the set-up ~1+condition+batch is perfectly confounded. + + Can only group by non-constrained effects right now, use constraint_matrix_from_string + for other cases. + - list of strings or tuple of strings: + String encoded equality constraints. + + E.g. ["batch1 + batch2 + batch3 = 0"] + - None: + No constraints are used, this is equivalent to using an identity matrix as a + constraint matrix. + + Note that np.ndarray encoded full constraint matrices are not supported here as the design + matrices are built within this function. + :param constraints_scale: Constraints for scale model. Can be following: + + - dict: + Every element of the dictionary corresponds to one set of equality constraints. + Each set has to be be an entry of the form {..., x: y, ...} + where x is the factor to be constrained and y is a factor by which levels of x are grouped + and then constrained. Set y="1" to constrain all levels of x to sum to one, + a single equality constraint. + + E.g.: {"batch": "condition"} Batch levels within each condition are constrained to sum to + zero. This is applicable if repeats of a an experiment within each condition + are independent so that the set-up ~1+condition+batch is perfectly confounded. + + Can only group by non-constrained effects right now, use constraint_matrix_from_string + for other cases. + - list of strings or tuple of strings: + String encoded equality constraints. + + E.g. ["batch1 + batch2 + batch3 = 0"] + - None: + No constraints are used, this is equivalent to using an identity matrix as a + constraint matrix. + + Note that np.ndarray encoded full constraint matrices are not supported here as the design + matrices are built within this function. :param noise_model: str, noise model to use in model-based unit_test. Possible options: - 'nb': default @@ -1732,9 +1869,9 @@ def continuous_1d( - str: will use Estimator.TrainingStrategy[training_strategy] to train - function: Can be used to implement custom training function will be called as - `training_strategy(estimator)`. - - list of keyword dicts containing method arguments: Will call Estimator.train() once with each dict of - method arguments. + `training_strategy(estimator)`. + - list of keyword dicts containing method arguments: Will call Estimator.train() + once with each dict of method arguments. Example: @@ -1754,7 +1891,7 @@ def continuous_1d( Should be "float32" for single precision or "float64" for double precision. :param kwargs: [Debugging] Additional arguments will be passed to the _fit method. """ - if formula_loc is None: + if formula_loc is None: raise ValueError("supply fomula_loc") # Set testing default to continuous covariate if not supplied: if factor_loc_totest is None: @@ -1824,10 +1961,10 @@ def continuous_1d( else: factor_loc_totest_new = factor_loc_totest - logger.debug("model formulas assembled in de.test.continuos():") - logger.debug("factor_loc_totest_new: " + ",".join(factor_loc_totest_new)) - logger.debug("formula_loc_new: " + formula_loc_new) - logger.debug("formula_scale_new: " + formula_scale_new) + logging.getLogger("diffxpy").debug("model formulas assembled in de.test.continuos():") + logging.getLogger("diffxpy").debug("factor_loc_totest_new: " + ",".join(factor_loc_totest_new)) + logging.getLogger("diffxpy").debug("formula_loc_new: " + formula_loc_new) + logging.getLogger("diffxpy").debug("formula_scale_new: " + formula_scale_new) de_test = wald( data=X, @@ -1840,6 +1977,8 @@ def continuous_1d( init_b=init_b, gene_names=gene_names, sample_description=sample_description, + constraints_loc=constraints_loc, + constraints_scale=constraints_scale, noise_model=noise_model, size_factors=size_factors, batch_size=batch_size, @@ -1850,6 +1989,7 @@ def continuous_1d( ) de_test = DifferentialExpressionTestWaldCont( de_test=de_test, + noise_model=noise_model, size_factors=size_factors, continuous_coords=sample_description[continuous].values, spline_coefs=new_coefs @@ -1873,11 +2013,11 @@ def continuous_1d( full_formula_scale = formula_scale_new reduced_formula_scale = formula_scale_new - logger.debug("model formulas assembled in de.test.continuous():") - logger.debug("full_formula_loc: " + full_formula_loc) - logger.debug("reduced_formula_loc: " + reduced_formula_loc) - logger.debug("full_formula_scale: " + full_formula_scale) - logger.debug("reduced_formula_scale: " + reduced_formula_scale) + logging.getLogger("diffxpy").debug("model formulas assembled in de.test.continuous():") + logging.getLogger("diffxpy").debug("full_formula_loc: " + full_formula_loc) + logging.getLogger("diffxpy").debug("reduced_formula_loc: " + reduced_formula_loc) + logging.getLogger("diffxpy").debug("full_formula_scale: " + full_formula_scale) + logging.getLogger("diffxpy").debug("reduced_formula_scale: " + reduced_formula_scale) de_test = lrt( data=X, @@ -1905,6 +2045,6 @@ def continuous_1d( spline_coefs=new_coefs ) else: - raise ValueError('base.continuous(): Parameter `test` not recognized.') + raise ValueError('continuous(): Parameter `test` not recognized.') - return de_test \ No newline at end of file + return de_test diff --git a/diffxpy/testing/utils.py b/diffxpy/testing/utils.py index 413a87a..3c348f1 100644 --- a/diffxpy/testing/utils.py +++ b/diffxpy/testing/utils.py @@ -1,17 +1,27 @@ -from typing import Union - import anndata import numpy as np import pandas as pd import patsy +import scipy +from typing import List, Tuple, Union import xarray as xr +try: + from anndata.base import Raw +except ImportError: + from anndata import Raw + from batchglm import data as data_utils +# Relay util functions for diffxpy api. +# design_matrix, preview_coef_names and constraint_system_from_star are redefined here. +from batchglm.data import constraint_matrix_from_string, constraint_matrix_from_dict +from batchglm.data import design_matrix_from_xarray, design_matrix_from_anndata +from batchglm.data import view_coef_names 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)): + if anndata is not None and (isinstance(data, anndata.AnnData) or isinstance(data, Raw)): gene_names = data.var_names elif isinstance(data, xr.DataArray): gene_names = data["features"] @@ -31,7 +41,17 @@ def parse_data(data, gene_names) -> xr.DataArray: return X -def parse_sample_description(data, sample_description=None) -> pd.DataFrame: +def parse_sample_description( + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + sample_description: Union[pd.DataFrame, None] +) -> pd.DataFrame: + """ + Parse sample description from input. + + :param data: Input data matrix (observations x features) or (cells x genes). + :param sample_description: pandas.DataFrame containing sample annotations, can be None. + :return: Assembled sample annotations. + """ if sample_description is None: if anndata is not None and isinstance(data, anndata.AnnData): sample_description = data_utils.sample_description_from_anndata( @@ -48,8 +68,8 @@ def parse_sample_description(data, sample_description=None) -> pd.DataFrame: "with corresponding sample annotations" ) - if anndata is not None and isinstance(data, anndata.base.Raw): - # anndata.base.Raw does not have attribute shape. + if anndata is not None and isinstance(data, Raw): + # 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: @@ -58,88 +78,207 @@ def parse_sample_description(data, sample_description=None) -> pd.DataFrame: return sample_description -def parse_size_factors(size_factors, data): +def parse_size_factors( + size_factors: Union[np.ndarray, pd.core.series.Series, np.ndarray], + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, scipy.sparse.csr_matrix], + sample_description: pd.DataFrame +) -> Union[np.ndarray, None]: + """ + Parse size-factors from input. + + :param size_factors: 1D array of transformed library size factors for each cell in the + same order as in data or string-type column identifier of size-factor containing + column in sample description. + :param data: Input data matrix (observations x features) or (cells x genes). + :param sample_description: optional pandas.DataFrame containing sample annotations + :return: Assebled size-factors. + """ if size_factors is not None: if isinstance(size_factors, pd.core.series.Series): size_factors = size_factors.values + elif isinstance(size_factors, str): + assert size_factors in sample_description.columns, "" + size_factors = sample_description[size_factors].values assert size_factors.shape[0] == data.shape[0], "data matrix and size factors must contain same number of cells" + assert np.all(size_factors > 0), "size_factors <= 0 found, please remove these cells" return size_factors +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 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 + + 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 + data: Union[anndata.AnnData, Raw, xr.DataArray, xr.Dataset, np.ndarray, + scipy.sparse.csr_matrix] = None, + sample_description: Union[None, pd.DataFrame] = None, + formula: Union[None, str] = None, + as_numeric: Union[List[str], Tuple[str], str] = (), + dmat: Union[pd.DataFrame, None] = None, + return_type: str = "xarray" +) -> Union[patsy.design_info.DesignMatrix, xr.Dataset, pd.DataFrame]: + """ Create a design matrix from some sample description. + + This function defaults to perform formatting if dmat is directly supplied as a pd.DataFrame. + This function relays batchglm.data.design_matrix() to behave like the other wrappers in diffxpy. + + :param data: Input data matrix (observations x features) or (cells x genes). + :param sample_description: pandas.DataFrame of length "num_observations" containing explanatory variables as columns + :param formula: model formula as string, describing the relations of the explanatory variables. + + E.g. '~ 1 + batch + confounder' + :param as_numeric: + Which columns of sample_description to treat as numeric and + not as categorical. This yields columns in the design matrix + which do not correpond to one-hot encoded discrete factors. + This makes sense for number of genes, time, pseudotime or space + for example. + :param dmat: a model design matrix as a pd.DataFrame + :param return_type: type of the returned value. + + - "patsy": return plain patsy.design_info.DesignMatrix object + - "dataframe": return pd.DataFrame with observations as rows and params as columns + - "xarray": return xr.Dataset with design matrix as ds["design"] and the sample description embedded as + one variable per column :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.") + 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.") + 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 + if sample_description is not None: + as_categorical = [False if x in as_numeric else True for x in sample_description.columns.values] else: - ar = xr.DataArray(dmat, dims=("observations", "design_params")) - ar.coords["design_params"] = dmat.columns + as_categorical = True + + return data_utils.design_matrix( + sample_description=sample_description, + formula=formula, + as_categorical=as_categorical, + dmat=dmat, + return_type=return_type + ) - ds = xr.Dataset({ - "design": ar, - }) - return ds +def preview_coef_names( + sample_description: pd.DataFrame, + formula: str, + as_numeric: Union[List[str], Tuple[str], str] = () +) -> np.ndarray: + """ + Return coefficient names of model. + Use this to preview what the model would look like. + This function relays batchglm.data.preview_coef_names() to behave like the other wrappers in diffxpy. -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 sample_description: pandas.DataFrame of length "num_observations" containing explanatory variables as columns + :param formula: model formula as string, describing the relations of the explanatory variables. - :param data: input data - :param formula: model formula. - :param sample_description: optional pandas.DataFrame containing sample annotations - :param dmat: model design matrix + E.g. '~ 1 + batch + confounder' + :param as_numeric: + Which columns of sample_description to treat as numeric and + not as categorical. This yields columns in the design matrix + which do not correpond to one-hot encoded discrete factors. + This makes sense for number of genes, time, pseudotime or space + for example. + :return: A list of coefficient names. """ - return design_matrix( - data=data, + if isinstance(as_numeric, str): + as_numeric = [as_numeric] + if isinstance(as_numeric, tuple): + as_numeric = list(as_numeric) + + return data_utils.preview_coef_names( sample_description=sample_description, formula=formula, - dmat=dmat - ).design_info.column_names + as_categorical=[False if x in as_numeric else True for x in sample_description.columns.values] + ) -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 constraint_system_from_star( + dmat: Union[None, np.ndarray, xr.DataArray, xr.Dataset] = None, + sample_description: Union[None, pd.DataFrame] = None, + formula: Union[None, str] = None, + as_numeric: Union[List[str], Tuple[str], str] = (), + constraints: dict = {}, + dims: Union[Tuple[str, str], List[str]] = (), + return_type: str = "xarray", +) -> Tuple: + """ + Create a design matrix and a constraint matrix. + This function relays batchglm.data.constraint_matrix_from_star() to behave like the other wrappers in diffxpy. -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 + :param dmat: Pre-built model design matrix. + :param sample_description: pandas.DataFrame of length "num_observations" containing explanatory variables as columns + :param formula: model formula as string, describing the relations of the explanatory variables. + E.g. '~ 1 + batch + confounder' + :param as_numeric: + Which columns of sample_description to treat as numeric and + not as categorical. This yields columns in the design matrix + which do not correspond to one-hot encoded discrete factors. + :param constraints: Grouped factors to enfore equality constraints on. Every element of + the dictionary corresponds to one set of equality constraints. Each set has to be + be an entry of the form {..., x: y, ...} where x is the factor to be constrained and y is + a factor by which levels of x are grouped and then constrained. Set y="1" to constrain + all levels of x to sum to one, a single equality constraint. -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) + E.g.: {"batch": "condition"} Batch levels within each condition are constrained to sum to + zero. This is applicable if repeats of a an experiment within each condition + are independent so that the set-up ~1+condition+batch is perfectly confounded. + + Can only group by non-constrained effects right now, use constraint_matrix_from_string + for other cases. + :param dims: Dimension names of xarray. - return dmat, sample_description \ No newline at end of file + E.g.: ["design_loc_params", "loc_params"] or ["design_scale_params", "scale_params"] + :param return_type: type of the returned value. + + - "patsy": return plain patsy.design_info.DesignMatrix object + - "dataframe": return pd.DataFrame with observations as rows and params as columns + - "xarray": return xr.Dataset with design matrix as ds["design"] and the sample description embedded as + one variable per column + This option is overridden if constraints are supplied as dict. + :return: a model design matrix and a constraint matrix formatted as xr.DataArray + """ + if isinstance(as_numeric, str): + as_numeric = [as_numeric] + if isinstance(as_numeric, tuple): + as_numeric = list(as_numeric) + + if sample_description is not None: + as_categorical = [False if x in as_numeric else True for x in sample_description.columns.values] + else: + as_categorical = True + + return data_utils.constraint_system_from_star( + dmat=dmat, + sample_description=sample_description, + formula=formula, + as_categorical=as_categorical, + constraints=constraints, + dims=dims, + return_type=return_type + ) diff --git a/diffxpy/unit_test/test_constrained.py b/diffxpy/unit_test/test_constrained.py index f674820..9ef21a3 100644 --- a/diffxpy/unit_test/test_constrained.py +++ b/diffxpy/unit_test/test_constrained.py @@ -9,24 +9,20 @@ import diffxpy.api as de -class TestSingle(unittest.TestCase): +class TestConstrained(unittest.TestCase): - def test_null_distribution_wald_constrained(self, n_genes: int = 100): + def test_forfatal_from_string(self): """ - Test if de.wald() with constraints generates a uniform p-value distribution - if it is given data simulated based on the null model. Returns the p-value - of the two-side Kolmgorov-Smirnov test for equality of the observed - p-value distribution and a uniform distribution. + Test if _from_string interface is working. n_cells is constant as the design matrix and constraints depend on it. - - :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) n_cells = 2000 + n_genes = 2 sim = Simulator(num_observations=n_cells, num_features=n_genes) sim.generate_sample_description(num_batches=0, num_conditions=0) @@ -43,16 +39,16 @@ def test_null_distribution_wald_constrained(self, n_genes: int = 100): coefficient_names = ['intercept', 'bio1', 'bio2', 'bio3', 'bio4', 'treatment1'] dmat_est = pd.DataFrame(data=dmat, columns=coefficient_names) - dmat_est_loc = de.test.design_matrix(dmat=dmat_est) - dmat_est_scale = de.test.design_matrix(dmat=dmat_est) + dmat_est_loc = de.utils.design_matrix(dmat=dmat_est) + dmat_est_scale = de.utils.design_matrix(dmat=dmat_est) # Build constraints: - constraints_loc = de.utils.data_utils.build_equality_constraints_string( + constraints_loc = de.utils.constraint_matrix_from_string( dmat=dmat_est_loc, constraints=["bio1+bio2=0", "bio3+bio4=0"], dims=["design_loc_params", "loc_params"] ) - constraints_scale = de.utils.data_utils.build_equality_constraints_string( + constraints_scale = de.utils.constraint_matrix_from_string( dmat=dmat_est_scale, constraints=["bio1+bio2=0", "bio3+bio4=0"], dims=["design_scale_params", "scale_params"] @@ -60,15 +56,111 @@ def test_null_distribution_wald_constrained(self, n_genes: int = 100): test = de.test.wald( data=sim.X, - dmat_loc=dmat_est_loc.data_vars['design'], - dmat_scale=dmat_est_scale.data_vars['design'], - init_a="standard", - init_b="standard", + dmat_loc=dmat_est_loc, + dmat_scale=dmat_est_scale, + constraints_loc=constraints_loc, + constraints_scale=constraints_scale, + coef_to_test=["treatment1"] + ) + summary = test.summary() + + return True + + def test_forfatal_from_dict(self): + """ + Test if dictionary-based constraint interface is working. + """ + logging.getLogger("tensorflow").setLevel(logging.ERROR) + logging.getLogger("batchglm").setLevel(logging.WARNING) + logging.getLogger("diffxpy").setLevel(logging.WARNING) + + n_cells = 2000 + n_genes = 2 + + sim = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=0) + sim.generate() + + # Build design matrix: + sample_description = pd.DataFrame({ + "cond": ["cond"+str(i // 1000) for i in range(n_cells)], + "batch": ["batch"+str(i // 500) for i in range(n_cells)] + }) + + # Build constraints: + dmat_loc, constraints_loc = de.utils.constraint_matrix_from_dict( + sample_description=sample_description, + formula="~1+cond+batch", + constraints={"batch": "cond"}, + dims=["design_loc_params", "loc_params"] + ) + dmat_scale, constraints_scale = de.utils.constraint_matrix_from_dict( + sample_description=sample_description, + formula="~1+cond+batch", + constraints={"batch": "cond"}, + dims=["design_scale_params", "scale_params"] + ) + + test = de.test.wald( + data=sim.X, + dmat_loc=dmat_loc, + dmat_scale=dmat_scale, constraints_loc=constraints_loc, constraints_scale=constraints_scale, - coef_to_test=["treatment1"], - training_strategy="DEFAULT", - dtype="float64" + coef_to_test=["cond[T.cond1]"] + ) + summary = test.summary() + + return True + + def test_null_distribution_wald_constrained(self, n_genes: int = 100): + """ + Test if de.wald() with constraints generates a uniform p-value distribution + if it is given data simulated based on the null model. Returns the p-value + of the two-side Kolmgorov-Smirnov test for equality of the observed + p-value distribution and a uniform distribution. + + n_cells is constant as the design matrix and constraints depend on it. + + :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) + + n_cells = 2000 + + sim = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=0) + sim.generate() + + # Build design matrix: + sample_description = pd.DataFrame({ + "cond": ["cond" + str(i // 1000) for i in range(n_cells)], + "batch": ["batch" + str(i // 500) for i in range(n_cells)] + }) + + # Build constraints: + dmat_loc, constraints_loc = de.utils.constraint_matrix_from_dict( + sample_description=sample_description, + formula="~1+cond+batch", + constraints={"batch": "cond"}, + dims=["design_loc_params", "loc_params"] + ) + dmat_scale, constraints_scale = de.utils.constraint_matrix_from_dict( + sample_description=sample_description, + formula="~1+cond+batch", + constraints={"batch": "cond"}, + dims=["design_scale_params", "scale_params"] + ) + + test = de.test.wald( + data=sim.X, + dmat_loc=dmat_loc, + dmat_scale=dmat_scale, + constraints_loc=constraints_loc, + constraints_scale=constraints_scale, + coef_to_test=["cond[T.cond1]"] ) summary = test.summary() @@ -127,11 +219,11 @@ def test_null_distribution_wald_constrained_2layer(self, n_genes: int = 100): 'tech1', 'tech2', 'tech3', 'tech4'] dmat_est = pd.DataFrame(data=dmat, columns=coefficient_names) - dmat_est_loc = de.test.design_matrix(dmat=dmat_est) - dmat_est_scale = de.test.design_matrix(dmat=dmat_est.iloc[:, [0]]) + dmat_est_loc = de.utils.design_matrix(dmat=dmat_est) + dmat_est_scale = de.utils.design_matrix(dmat=dmat_est.iloc[:, [0]]) # Build constraints: - constraints_loc = de.utils.data_utils.build_equality_constraints_string( + constraints_loc = de.utils.constraint_matrix_from_string( dmat=dmat_est_loc, constraints=["bio1+bio2=0", "bio3+bio4=0", @@ -145,16 +237,11 @@ def test_null_distribution_wald_constrained_2layer(self, n_genes: int = 100): test = de.test.wald( data=sim.X, - dmat_loc=dmat_est_loc.data_vars['design'], - dmat_scale=dmat_est_scale.data_vars['design'], - init_a="standard", - init_b="standard", + dmat_loc=dmat_est_loc, + dmat_scale=dmat_est_scale, constraints_loc=constraints_loc, constraints_scale=constraints_scale, - coef_to_test=["treatment1"], - training_strategy="DEFAULT", - quick_scale=False, - dtype="float64" + coef_to_test=["treatment1"] ) summary = test.summary() @@ -203,18 +290,18 @@ def test_null_distribution_wald_multi_constrained_2layer(self, n_genes: int = 50 'bio5', 'bio6', 'treatment1', 'treatment2'] dmat_est = pd.DataFrame(data=dmat, columns=coefficient_names) - dmat_est_loc = de.test.design_matrix(dmat=dmat_est) - dmat_est_scale = de.test.design_matrix(dmat=dmat_est) + dmat_est_loc = de.utils.design_matrix(dmat=dmat_est) + dmat_est_scale = de.utils.design_matrix(dmat=dmat_est) # Build constraints: - constraints_loc = de.utils.data_utils.build_equality_constraints_string( + constraints_loc = de.utils.constraint_matrix_from_string( dmat=dmat_est_loc, constraints=["bio1+bio2=0", "bio3+bio4=0", "bio5+bio6=0"], dims=["design_loc_params", "loc_params"] ) - constraints_scale = de.utils.data_utils.build_equality_constraints_string( + constraints_scale = de.utils.constraint_matrix_from_string( dmat=dmat_est_scale, constraints=["bio1+bio2=0", "bio3+bio4=0", @@ -224,13 +311,11 @@ def test_null_distribution_wald_multi_constrained_2layer(self, n_genes: int = 50 test = de.test.wald( data=sim.X, - dmat_loc=dmat_est_loc.data_vars['design'], - dmat_scale=dmat_est_scale.data_vars['design'], + dmat_loc=dmat_est_loc, + dmat_scale=dmat_est_scale, constraints_loc=constraints_loc, constraints_scale=constraints_scale, - coef_to_test=["treatment1", "treatment2"], - training_strategy="DEFAULT", - dtype="float64" + coef_to_test=["treatment1", "treatment2"] ) summary = test.summary() @@ -238,7 +323,7 @@ def test_null_distribution_wald_multi_constrained_2layer(self, n_genes: int = 50 pval_h0 = stats.kstest(test.pval, 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % pval_h0 return True diff --git a/diffxpy/unit_test/test_enrich.py b/diffxpy/unit_test/test_enrich.py new file mode 100644 index 0000000..8ed229f --- /dev/null +++ b/diffxpy/unit_test/test_enrich.py @@ -0,0 +1,57 @@ +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 TestEnrich(unittest.TestCase): + + def test_for_fatal(self): + """ + """ + logging.getLogger("tensorflow").setLevel(logging.ERROR) + logging.getLogger("batchglm").setLevel(logging.WARNING) + logging.getLogger("diffxpy").setLevel(logging.WARNING) + + sim = Simulator(num_observations=50, num_features=10) + sim.generate_sample_description(num_batches=0, num_conditions=2) + sim.generate() + + test = de.test.wald( + data=sim.X, + factor_loc_totest="condition", + formula_loc="~ 1 + condition", + sample_description=sim.sample_description, + gene_names=[str(x) for x in range(sim.X.shape[1])], + training_strategy="DEFAULT", + dtype="float64" + ) + + # Set up reference gene sets. + rs = de.enrich.RefSets() + rs.add(id="set1", source="manual", gene_ids=["1", "3"]) + rs.add(id="set2", source="manual", gene_ids=["5", "6"]) + + for i in [True, False]: + for j in [True, False]: + enrich_test_i = de.enrich.test( + ref=rs, + det=test, + threshold=0.05, + incl_all_zero=i, + clean_ref=j, + ) + _ = enrich_test_i.summary() + _ = enrich_test_i.significant_set_ids() + _ = enrich_test_i.significant_sets() + _ = enrich_test_i.set_summary(id="set1") + + return True + + +if __name__ == '__main__': + unittest.main() diff --git a/diffxpy/unit_test/test_extreme_values.py b/diffxpy/unit_test/test_extreme_values.py index f04d535..cb19922 100644 --- a/diffxpy/unit_test/test_extreme_values.py +++ b/diffxpy/unit_test/test_extreme_values.py @@ -11,24 +11,19 @@ class TestExtremeValues(unittest.TestCase): - def test_t_test_zero_variance(self, n_cells: int = 2000, n_genes: int = 100): + def test_t_test_zero_variance(self): """ - Test if de.t_test() generates a uniform p-value distribution - if it is given data simulated based on the null model. Returns the p-value - of the two-side Kolmgorov-Smirnov test for equality of the observed - p-value distribution and a uniform distribution. - - :param n_cells: Number of cells to simulate (number of observations per test). - :param n_genes: Number of genes to simulate (number of tests). + Test if T-test works if it is given genes with zero variance. """ logging.getLogger("tensorflow").setLevel(logging.ERROR) logging.getLogger("batchglm").setLevel(logging.WARNING) logging.getLogger("diffxpy").setLevel(logging.WARNING) - sim = Simulator(num_observations=n_cells, num_features=n_genes) + sim = Simulator(num_observations=1000, num_features=10) sim.generate_sample_description(num_batches=0, num_conditions=0) sim.generate() - sim.data.X[:, 0] = np.exp(sim.a)[0, 0] + sim.data.X[:, 0] = 0 + sim.data.X[:, 1] = 5 random_sample_description = pd.DataFrame({ "condition": np.random.randint(2, size=sim.num_observations) @@ -37,17 +32,44 @@ def test_t_test_zero_variance(self, n_cells: int = 2000, n_genes: int = 100): test = de.test.t_test( data=sim.X, grouping="condition", - sample_description=random_sample_description + sample_description=random_sample_description, + is_sig_zerovar=True ) - # Compare p-value distribution under null model against uniform distribution. - pval_h0 = stats.kstest(test.pval, 'uniform').pvalue + assert np.isnan(test.pval[0]) and test.pval[1] == 1, \ + "rank test did not assign p-value of zero to groups with zero variance and same mean, %f, %f" % \ + (test.pval[0], test.pval[1]) + return True + + def test_rank_test_zero_variance(self): + """ + Test if rank test works if it is given genes with zero variance. + """ + logging.getLogger("tensorflow").setLevel(logging.ERROR) + logging.getLogger("batchglm").setLevel(logging.WARNING) + logging.getLogger("diffxpy").setLevel(logging.WARNING) + + sim = Simulator(num_observations=1000, num_features=10) + sim.generate_sample_description(num_batches=0, num_conditions=0) + sim.generate() + sim.data.X[:, 0] = 0 + sim.data.X[:, 1] = 5 - print('KS-test pvalue for null model match of t_test(): %f' % pval_h0) + random_sample_description = pd.DataFrame({ + "condition": np.random.randint(2, size=sim.num_observations) + }) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + test = de.test.rank_test( + data=sim.X, + grouping="condition", + sample_description=random_sample_description, + is_sig_zerovar=True + ) - return pval_h0 + assert np.isnan(test.pval[0]) and test.pval[1] == 1, \ + "rank test did not assign p-value of zero to groups with zero variance and same mean, %f, %f" % \ + (test.pval[0], test.pval[1]) + return True if __name__ == '__main__': diff --git a/diffxpy/unit_test/test_pairwise.py b/diffxpy/unit_test/test_pairwise.py index 5e373fb..5b52746 100644 --- a/diffxpy/unit_test/test_pairwise.py +++ b/diffxpy/unit_test/test_pairwise.py @@ -46,7 +46,7 @@ def test_null_distribution_ztest(self, n_cells: int = 2000, n_genes: int = 100, pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True @@ -89,7 +89,7 @@ def test_null_distribution_z_lazy(self, n_cells: int = 2000, n_genes: int = 100) pval_h0 = stats.kstest(pvals.flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True @@ -128,7 +128,7 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100, n_ pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True @@ -166,7 +166,7 @@ def test_null_distribution_ttest(self, n_cells: int = 2000, n_genes: int = 10000 pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True @@ -204,7 +204,7 @@ def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 10 pval_h0 = stats.kstest(test.pval[~np.eye(test.pval.shape[0]).astype(bool)].flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True diff --git a/diffxpy/unit_test/test_partition.py b/diffxpy/unit_test/test_partition.py new file mode 100644 index 0000000..98d1009 --- /dev/null +++ b/diffxpy/unit_test/test_partition.py @@ -0,0 +1,237 @@ +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 TestPartitionNull(unittest.TestCase): + + def test_null_distribution_wald(self, n_cells: int = 4000, n_genes: int = 200): + """ + Test if Partition.wald() generates a uniform p-value distribution + if it is given data simulated based on the null model. Returns the p-value + of the two-side Kolmgorov-Smirnov test for equality of the observed + p-value distribution and a uniform distribution. + + :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 = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=2) + sim.generate() + + sample_description = pd.DataFrame({ + "covar1": np.random.randint(2, size=sim.num_observations), + "covar2": np.random.randint(2, size=sim.num_observations) + }) + sample_description["cond"] = sim.sample_description["condition"].values + + partition = de.test.partition( + data=sim.X, + parts="cond", + sample_description=sample_description + ) + det = partition.wald( + factor_loc_totest="covar1", + formula_loc="~ 1 + covar1 + covar2", + training_strategy="DEFAULT", + dtype="float64" + ) + summary = det.summary() + + # Compare p-value distribution under null model against uniform distribution. + pval_h0 = stats.kstest(det.pval.flatten(), 'uniform').pvalue + + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + + return True + + def test_null_distribution_wald_multi(self, n_cells: int = 4000, n_genes: int = 200): + """ + Test if de.wald() (multivariate mode) generates a uniform p-value distribution + if it is given data simulated based on the null model. Returns the p-value + of the two-side Kolmgorov-Smirnov test for equality of the observed + p-value distribution and a uniform distribution. + + :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 = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=2) + sim.generate() + + sample_description = pd.DataFrame({ + "covar1": np.random.randint(4, size=sim.num_observations), + "covar2": np.random.randint(2, size=sim.num_observations) + }) + sample_description["cond"] = sim.sample_description["condition"].values + + partition = de.test.partition( + data=sim.X, + parts="cond", + sample_description=sample_description + ) + det = partition.wald( + factor_loc_totest="covar1", + formula_loc="~ 1 + covar1 + covar2", + training_strategy="DEFAULT", + dtype="float64" + ) + summary = det.summary() + + # Compare p-value distribution under null model against uniform distribution. + pval_h0 = stats.kstest(det.pval.flatten(), 'uniform').pvalue + + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + + return True + + def test_null_distribution_lrt(self, n_cells: int = 4000, n_genes: int = 200): + """ + Test if de.lrt() generates a uniform p-value distribution + if it is given data simulated based on the null model. Returns the p-value + of the two-side Kolmgorov-Smirnov test for equality of the observed + p-value distribution and a uniform distribution. + + :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 = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=2) + sim.generate() + + sample_description = pd.DataFrame({ + "covar1": np.random.randint(2, size=sim.num_observations), + "covar2": np.random.randint(2, size=sim.num_observations) + }) + sample_description["cond"] = sim.sample_description["condition"].values + + partition = de.test.partition( + data=sim.X, + parts="cond", + sample_description=sample_description + ) + det = partition.lrt( + full_formula_loc="~ 1 + covar1", + full_formula_scale="~ 1", + reduced_formula_loc="~ 1", + reduced_formula_scale="~ 1", + training_strategy="DEFAULT", + dtype="float64" + ) + summary = det.summary() + + # Compare p-value distribution under null model against uniform distribution. + pval_h0 = stats.kstest(det.pval.flatten(), 'uniform').pvalue + + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of lrt(): %f' % pval_h0) + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + + return True + + def test_null_distribution_ttest(self, n_cells: int = 4000, n_genes: int = 200): + """ + Test if de.t_test() generates a uniform p-value distribution + if it is given data simulated based on the null model. Returns the p-value + of the two-side Kolmgorov-Smirnov test for equality of the observed + p-value distribution and a uniform distribution. + + :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 = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=2) + sim.generate() + + sample_description = pd.DataFrame({ + "covar1": np.random.randint(2, size=sim.num_observations) + }) + sample_description["cond"] = sim.sample_description["condition"].values + + partition = de.test.partition( + data=sim.X, + parts="cond", + sample_description=sample_description + ) + det = partition.t_test( + grouping="covar1", + is_logged=False, + dtype="float64" + ) + summary = det.summary() + + # Compare p-value distribution under null model against uniform distribution. + pval_h0 = stats.kstest(det.pval.flatten(), 'uniform').pvalue + + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of t_test(): %f' % pval_h0) + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + + return True + + def test_null_distribution_rank(self, n_cells: int = 4000, n_genes: int = 200): + """ + Test if rank_test() generates a uniform p-value distribution + if it is given data simulated based on the null model. Returns the p-value + of the two-side Kolmgorov-Smirnov test for equality of the observed + p-value distribution and a uniform distribution. + + :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 = Simulator(num_observations=n_cells, num_features=n_genes) + sim.generate_sample_description(num_batches=0, num_conditions=2) + sim.generate() + + sample_description = pd.DataFrame({ + "covar1": np.random.randint(2, size=sim.num_observations) + }) + sample_description["cond"] = sim.sample_description["condition"].values + + partition = de.test.partition( + data=sim.X, + parts="cond", + sample_description=sample_description + ) + det = partition.rank_test( + grouping="covar1", + dtype="float64" + ) + summary = det.summary() + + # Compare p-value distribution under null model against uniform distribution. + pval_h0 = stats.kstest(det.pval.flatten(), 'uniform').pvalue + + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of rank_test(): %f' % pval_h0) + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + + return True + + +if __name__ == '__main__': + unittest.main() diff --git a/diffxpy/unit_test/test_single_de.py b/diffxpy/unit_test/test_single_de.py new file mode 100644 index 0000000..693e1fc --- /dev/null +++ b/diffxpy/unit_test/test_single_de.py @@ -0,0 +1,305 @@ +import unittest +import logging +import numpy as np +import pandas as pd +import scipy.stats as stats + +import diffxpy.api as de + + +class _TestSingleDE: + + def _prepare_data( + 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. + """ + if noise_model == "nb": + from batchglm.api.models.glm_nb import Simulator + elif noise_model == "norm": + from batchglm.api.models.glm_norm import Simulator + else: + raise ValueError("noise model %s not recognized" % noise_model) + + 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, + n_genes: int + ): + """ + :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, + noise_model="norm" + ) + + 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, + n_genes: int + ): + """ + :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, + noise_model="norm" + ) + + 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, + noise_model=noise_model + ) + + 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, + noise_model=noise_model + ) + + 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 new file mode 100644 index 0000000..0590c2a --- /dev/null +++ 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.py b/diffxpy/unit_test/test_single_null.py similarity index 53% rename from diffxpy/unit_test/test_single.py rename to diffxpy/unit_test/test_single_null.py index 37e760f..34ca47d 100644 --- a/diffxpy/unit_test/test_single.py +++ b/diffxpy/unit_test/test_single_null.py @@ -4,13 +4,17 @@ import pandas as pd import scipy.stats as stats -from batchglm.api.models.glm_nb import Simulator import diffxpy.api as de -class TestSingleNull(unittest.TestCase): +class _TestSingleNull: - def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 100): + def _test_null_distribution_wald( + self, + n_cells: int, + n_genes: int, + noise_model: str + ): """ Test if de.wald() generates a uniform p-value distribution if it is given data simulated based on the null model. Returns the p-value @@ -19,10 +23,14 @@ def test_null_distribution_wald(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). + :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) + if noise_model == "nb": + from batchglm.api.models.glm_nb import Simulator + elif noise_model == "norm": + from batchglm.api.models.glm_norm import Simulator + else: + raise ValueError("noise model %s not recognized" % noise_model) sim = Simulator(num_observations=n_cells, num_features=n_genes) sim.generate_sample_description(num_batches=0, num_conditions=0) @@ -39,6 +47,7 @@ def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 100): formula_loc="~ 1 + condition + batch", sample_description=random_sample_description, batch_size=500, + noise_model=noise_model, training_strategy="DEFAULT", dtype="float64" ) @@ -48,11 +57,16 @@ def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 100): pval_h0 = stats.kstest(test.pval, 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + assert pval_h0 > 0.05, ("KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)) return True - def test_null_distribution_wald_multi(self, n_cells: int = 2000, n_genes: int = 100): + def _test_null_distribution_wald_multi( + self, + n_cells: int, + n_genes: int, + noise_model: str + ): """ Test if de.wald() (multivariate mode) generates a uniform p-value distribution if it is given data simulated based on the null model. Returns the p-value @@ -61,10 +75,14 @@ def test_null_distribution_wald_multi(self, n_cells: int = 2000, n_genes: int = :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) + if noise_model == "nb": + from batchglm.api.models.glm_nb import Simulator + elif noise_model == "norm": + from batchglm.api.models.glm_norm import Simulator + else: + raise ValueError("noise model %s not recognized" % noise_model) sim = Simulator(num_observations=n_cells, num_features=n_genes) sim.generate_sample_description(num_batches=0, num_conditions=0) @@ -79,6 +97,7 @@ def test_null_distribution_wald_multi(self, n_cells: int = 2000, n_genes: int = factor_loc_totest="condition", formula_loc="~ 1 + condition", sample_description=random_sample_description, + noise_model=noise_model, training_strategy="DEFAULT", dtype="float64" ) @@ -88,11 +107,16 @@ def test_null_distribution_wald_multi(self, n_cells: int = 2000, n_genes: int = pval_h0 = stats.kstest(test.pval, 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + assert pval_h0 > 0.05, ("KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)) return True - def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100): + def _test_null_distribution_lrt( + self, + n_cells: int, + n_genes: int, + noise_model: str + ): """ Test if de.lrt() generates a uniform p-value distribution if it is given data simulated based on the null model. Returns the p-value @@ -101,10 +125,14 @@ def test_null_distribution_lrt(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). + :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) + if noise_model == "nb": + from batchglm.api.models.glm_nb import Simulator + elif noise_model == "norm": + from batchglm.api.models.glm_norm import Simulator + else: + raise ValueError("noise model %s not recognized" % noise_model) sim = Simulator(num_observations=n_cells, num_features=n_genes) sim.generate_sample_description(num_batches=0, num_conditions=0) @@ -121,6 +149,7 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100): reduced_formula_loc="~ 1", reduced_formula_scale="~ 1", sample_description=random_sample_description, + noise_model=noise_model, training_strategy="DEFAULT", dtype="float64" ) @@ -130,11 +159,15 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100): pval_h0 = stats.kstest(test.pval, 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of lrt(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + assert pval_h0 > 0.05, ("KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)) return True - def test_null_distribution_ttest(self, n_cells: int = 2000, n_genes: int = 100): + def _test_null_distribution_ttest( + self, + n_cells: int, + n_genes: int + ): """ Test if de.t_test() generates a uniform p-value distribution if it is given data simulated based on the null model. Returns the p-value @@ -144,9 +177,7 @@ def test_null_distribution_ttest(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) + from batchglm.api.models.glm_norm import Simulator sim = Simulator(num_observations=n_cells, num_features=n_genes) sim.generate_sample_description(num_batches=0, num_conditions=0) @@ -169,13 +200,17 @@ def test_null_distribution_ttest(self, n_cells: int = 2000, n_genes: int = 100): pval_h0 = stats.kstest(test.pval, 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of t_test(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + assert pval_h0 > 0.05, ("KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)) return True - def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 100): + def _test_null_distribution_rank( + self, + n_cells: int, + n_genes: int + ): """ - Test if de.wilcoxon() generates a uniform p-value distribution + Test if de.test.rank_test() generates a uniform p-value distribution if it is given data simulated based on the null model. Returns the p-value of the two-side Kolmgorov-Smirnov test for equality of the observed p-value distribution and a uniform distribution. @@ -183,9 +218,7 @@ def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 10 :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) + from batchglm.api.models.glm_norm import Simulator sim = Simulator(num_observations=n_cells, num_features=n_genes) sim.generate_sample_description(num_batches=0, num_conditions=0) @@ -206,58 +239,44 @@ def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 10 # Compare p-value distribution under null model against uniform distribution. pval_h0 = stats.kstest(test.pval, 'uniform').pvalue - logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wilcoxon(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of rank_test(): %f' % pval_h0) + assert pval_h0 > 0.05, ("KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5)) return True -class TestSingleDE(unittest.TestCase): - - def _prepare_data(self, n_cells: int = 2000, n_genes: int = 100): +class TestSingleNullStandard(_TestSingleNull, unittest.TestCase): + """ + Noise model-independent unit tests that test whether a test generates uniformly + distributed p-values if data are sampled from the null model. + """ + def test_null_distribution_ttest( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ + Test if t_test() generates a uniform p-value distribution. :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("tensorflow").setLevel(logging.ERROR) + logging.getLogger("batchglm").setLevel(logging.WARNING) + logging.getLogger("diffxpy").setLevel(logging.WARNING) - logging.getLogger("diffxpy").info( - 'fraction of non-DE genes with q-value < 0.05: %.1f%%' % - float(100 * frac_de_of_non_de) + return self._test_null_distribution_ttest( + n_cells=n_cells, + n_genes=n_genes ) - 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_wilcoxon_de(self, n_cells: int = 2000, n_genes: int = 100): + def test_null_distribution_rank( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ - Test if de.test.t_test() generates a uniform p-value distribution - if it is given data simulated based on the null model. + Test if t_test() generates a uniform p-value distribution. :param n_cells: Number of cells to simulate (number of observations per test). :param n_genes: Number of genes to simulate (number of tests). @@ -266,23 +285,25 @@ def test_wilcoxon_de(self, n_cells: int = 2000, n_genes: int = 100): 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" + return self._test_null_distribution_rank( + n_cells=n_cells, + n_genes=n_genes ) - self._eval(sim=sim, test=test) - return True +class TestSingleNullNB(_TestSingleNull, unittest.TestCase): + """ + Negative binomial noise model unit tests that test whether a test generates uniformly + distributed p-values if data are sampled from the null model. + """ - def test_t_test_de(self, n_cells: int = 2000, n_genes: int = 100): + def test_null_distribution_wald_nb( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ - Test if de.test.t_test() generates a uniform p-value distribution - if it is given data simulated based on the null model. + Test if wald() generates a uniform p-value distribution for "nb" noise model. :param n_cells: Number of cells to simulate (number of observations per test). :param n_genes: Number of genes to simulate (number of tests). @@ -291,23 +312,20 @@ def test_t_test_de(self, n_cells: int = 2000, n_genes: int = 100): 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" + return self._test_null_distribution_wald( + n_cells=n_cells, + n_genes=n_genes, + noise_model="nb" ) - self._eval(sim=sim, test=test) - - return True - - def test_wald_de(self, n_cells: int = 2000, n_genes: int = 100): + def test_null_distribution_wald_multi_nb( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ - Test if de.test.wald() generates a uniform p-value distribution - if it is given data simulated based on the null model. + Test if wald() generates a uniform p-value distribution for "nb" noise model + for multiple coefficients to 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). @@ -316,27 +334,19 @@ def test_wald_de(self, n_cells: int = 2000, n_genes: int = 100): 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, - training_strategy="DEFAULT", - dtype="float64" + return self._test_null_distribution_wald_multi( + n_cells=n_cells, + n_genes=n_genes, + noise_model="nb" ) - self._eval(sim=sim, test=test) - - return True - - def test_lrt_de(self, n_cells: int = 2000, n_genes: int = 100): + def test_null_distribution_lrt_nb( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ - Test if de.test.lrt() generates a uniform p-value distribution - if it is given data simulated based on the null model. Returns the p-value - of the two-side Kolmgorov-Smirnov test for equality of the observed - p-value distribution and a uniform distribution. + Test if lrt() generates a uniform p-value distribution for "nb" noise model. :param n_cells: Number of cells to simulate (number of observations per test). :param n_genes: Number of genes to simulate (number of tests). @@ -345,128 +355,81 @@ def test_lrt_de(self, n_cells: int = 2000, n_genes: int = 100): 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, - training_strategy="DEFAULT", - dtype="float64" + return self._test_null_distribution_lrt( + n_cells=n_cells, + n_genes=n_genes, + noise_model="nb" ) - self._eval(sim=sim, test=test) - - return True - - -class TestSingleExternal(unittest.TestCase): - def _prepare_data(self, n_cells: int = 2000, n_genes: int = 100): +class TestSingleNullNORM(_TestSingleNull, unittest.TestCase): + """ + Normal noise model unit tests that test whether a test generates uniformly + distributed p-values if data are sampled from the null model. + """ + def test_null_distribution_wald_norm( + self, + n_cells: int = 200, + n_genes: int = 200 + ): """ + Test if wald() generates a uniform p-value distribution for "norm" noise model. :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) + logging.getLogger("tensorflow").setLevel(logging.ERROR) + logging.getLogger("batchglm").setLevel(logging.WARNING) + logging.getLogger("diffxpy").setLevel(logging.WARNING) + + return self._test_null_distribution_wald( + n_cells=n_cells, + n_genes=n_genes, + noise_model="norm" ) - 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): + def test_null_distribution_wald_multi_norm( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ - Test if de.test.t_test() generates the same p-value distribution as scipy t-test. + Test if wald() generates a uniform p-value distribution for "norm" noise model + for multiple coefficients to 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) + logging.getLogger("diffxpy").setLevel(logging.WARNING) - test = de.test.t_test( - data=sim.X, - grouping="condition", - sample_description=sim.sample_description, - dtype="float64" + return self._test_null_distribution_wald_multi( + n_cells=n_cells, + n_genes=n_genes, + noise_model="norm" ) - # 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_wilcoxon_ref(self, n_cells: int = 2000, n_genes: int = 100): + def test_null_distribution_lrt_norm( + self, + n_cells: int = 2000, + n_genes: int = 200 + ): """ - Test if de.test.t_test() generates the same p-value distribution as scipy t-test. + Test if lrt() generates a uniform p-value distribution for "norm" noise model. :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) + logging.getLogger("diffxpy").setLevel(logging.WARNING) - test = de.test.rank_test( - data=sim.X, - grouping="condition", - sample_description=sim.sample_description, - dtype="float64" + return self._test_null_distribution_lrt( + n_cells=n_cells, + n_genes=n_genes, + noise_model="norm" ) - # 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_vsrest.py b/diffxpy/unit_test/test_vsrest.py index 594fb51..28a56cf 100644 --- a/diffxpy/unit_test/test_vsrest.py +++ b/diffxpy/unit_test/test_vsrest.py @@ -48,7 +48,7 @@ def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 100, n pval_h0 = stats.kstest(test.pval.flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of test_wald_loc(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True @@ -64,7 +64,7 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100): """ logging.getLogger("tensorflow").setLevel(logging.ERROR) logging.getLogger("batchglm").setLevel(logging.WARNING) - logging.getLogger("diffxpy").setLevel(logging.WARNING) + logging.getLogger("diffxpy").setLevel(logging.ERROR) from batchglm.api.models.glm_nb import Simulator sim = Simulator(num_observations=n_cells, num_features=n_genes) @@ -90,12 +90,12 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100): # Compare p-value distribution under null model against uniform distribution. pval_h0 = stats.kstest(test.pval.flatten(), 'uniform').pvalue - logging.getLogger("diffxpy").info('KS-test pvalue for null model match of test_wald_loc(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + logging.getLogger("diffxpy").info('KS-test pvalue for null model match of test_lrt(): %f' % pval_h0) + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True - def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 100, n_groups: int = 2): + def test_null_distribution_rank(self, n_cells: int = 2000, n_genes: int = 100, n_groups: int = 2): """ Test if de.test_wald_loc() generates a uniform p-value distribution if it is given data simulated based on the null model. Returns the p-value @@ -121,7 +121,7 @@ def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 10 test = de.test.versus_rest( data=sim.X, grouping="condition", - test="wilcoxon", + test="rank", sample_description=random_sample_description, dtype="float64" ) @@ -131,7 +131,7 @@ def test_null_distribution_wilcoxon(self, n_cells: int = 2000, n_genes: int = 10 pval_h0 = stats.kstest(test.pval.flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of test_wald_loc(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True @@ -171,7 +171,7 @@ def test_null_distribution_ttest(self, n_cells: int = 2000, n_genes: int = 100, pval_h0 = stats.kstest(test.pval.flatten(), 'uniform').pvalue logging.getLogger("diffxpy").info('KS-test pvalue for null model match of test_wald_loc(): %f' % pval_h0) - assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!" + assert pval_h0 > 0.05, "KS-Test failed: pval_h0=%f is <= 0.05!" % np.round(pval_h0, 5) return True diff --git a/docs/api/index.rst b/docs/api/index.rst index 44e77cf..9a6345f 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -27,7 +27,7 @@ diffxpy provies infrastructure for likelihood ratio tests, Wald tests, t-tests a test.wald test.lrt test.t_test - test.wilcoxon + test.rank_test Multiple tests per gene ~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/conf.py b/docs/conf.py index 1e88cfc..09d9263 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,332 +1,161 @@ +# This code was adapted from https://github.com/theislab/scanpy/scanpy/conf.py +# This file is therefore licensed under the license of the scanpy project, +# available from https://github.com/theislab/scanpy and copied here at the time of accession. +# Note that multiple changes were made to this file to adapt it to the diffxpy project. + +# BSD 3-Clause License +# +# Copyright (c) 2017 F. Alexander Wolf, P. Angerer, Theis Lab +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# * Neither the name of the copyright holder nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + import sys -import inspect -import logging from pathlib import Path from datetime import datetime -from typing import Optional - -from sphinx.application import Sphinx -from sphinx.ext import autosummary - -# remove PyCharm’s old six module -if 'six' in sys.modules: - print(*sys.path, sep='\n') - for pypath in list(sys.path): - if any(p in pypath for p in ['PyCharm', 'pycharm']) and 'helpers' in pypath: - sys.path.remove(pypath) - del sys.modules['six'] - -import matplotlib # noqa -# Don’t use tkinter agg when importing scanpy → … → matplotlib +import matplotlib matplotlib.use('agg') HERE = Path(__file__).parent sys.path.insert(0, str(HERE.parent)) -import diffxpy # noqa +import diffxpy -logger = logging.getLogger(__name__) # -- General configuration ------------------------------------------------ -needs_sphinx = '1.7' # autosummary bugfix + +needs_sphinx = '1.7' + +# General information +project = 'diffxpy' +author = diffxpy.__author__ +copyright = f'{datetime.now():%Y}, {author}.' +version = diffxpy.__version__.replace('.dirty', '') +release = version + +# default settings +templates_path = ['_templates'] +source_suffix = '.rst' +master_doc = 'index' +default_role = 'literal' +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] +pygments_style = 'sphinx' + extensions = [ 'sphinx.ext.autodoc', + 'sphinx.ext.intersphinx', 'sphinx.ext.doctest', 'sphinx.ext.coverage', 'sphinx.ext.mathjax', - 'sphinx.ext.autosummary', - # 'plot_generator', - # 'plot_directive', 'sphinx.ext.napoleon', - 'sphinx_autodoc_typehints', - 'sphinx.ext.intersphinx', - # 'ipython_directive', - # 'ipython_console_highlighting', + 'sphinx.ext.autosummary' ] # Generate the API documentation when building autosummary_generate = True -# both of the following two lines don't work -# see falexwolf's issue for numpydoc -# autodoc_member_order = 'bysource' -# autodoc_default_flags = ['members'] +autodoc_member_order = 'bysource' napoleon_google_docstring = False napoleon_numpy_docstring = True napoleon_include_init_with_doc = False +napoleon_use_rtype = True +napoleon_use_param = True +napoleon_custom_sections = [('Params', 'Parameters')] +todo_include_todos = False intersphinx_mapping = dict( - python=('https://docs.python.org/3', None), + anndata=('https://anndata.readthedocs.io/en/latest/', None), + scanpy=('https://scanpy.readthedocs.io/en/latest/', None), numpy=('https://docs.scipy.org/doc/numpy/', None), - scipy=('https://docs.scipy.org/doc/scipy/reference/', None), pandas=('http://pandas.pydata.org/pandas-docs/stable/', None), - matplotlib=('https://matplotlib.org/', None), - # anndata=('https://anndata.readthedocs.io/en/latest/', None), + python=('https://docs.python.org/3', None), + scipy=('https://docs.scipy.org/doc/scipy/reference/', None) ) -templates_path = ['_templates'] - -project = 'diffxpy' -author = 'Florian R. Hölzlwimmer, David S. Fischer' - -source_suffix = '.rst' -master_doc = 'index' -copyright = f'{datetime.now():%Y}, {author}' - -version = diffxpy.__version__.replace('.dirty', '') -release = version -exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] -pygments_style = 'sphinx' -todo_include_todos = False # -- Options for HTML output ---------------------------------------------- + html_theme = 'sphinx_rtd_theme' html_theme_options = dict( - navigation_depth=2, + navigation_depth=4, + logo_only=True, # Only show the logo ) html_context = dict( - display_github=True, # Integrate GitHub - github_user='theislab', # Username - github_repo='diffxpy', # Repo name + display_github=True, # Integrate GitHub + github_user='theislab', # Username + github_repo='diffxpy', # Repo name github_version='master', # Version - conf_py_path='/docs/', # Path in the checkout to the docs root + conf_py_path='/docs/', # Path in the checkout to the docs root ) html_static_path = ['_static'] +html_show_sphinx = False +gh_url = '/{github_user}/{github_repo}'.format_map(html_context) def setup(app): app.add_stylesheet('css/custom.css') + app.connect('autodoc-process-docstring', insert_function_images) + app.add_role('pr', autolink(f'{gh_url}/pull/{{}}', 'PR {}')) -# -- Options for HTMLHelp output --------------------------------------------- - -# Output file base name for HTML help builder. -htmlhelp_basename = 'diffxpydoc' +# -- Options for other output formats ------------------------------------------ -# -- Options for LaTeX output ------------------------------------------------ -latex_elements = { - # The paper size ('letterpaper' or 'a4paper'). - # - # 'papersize': 'letterpaper', - - # The font size ('10pt', '11pt' or '12pt'). - # - # 'pointsize': '10pt', - - # Additional stuff for the LaTeX preamble. - # - # 'preamble': '', - - # Latex figure (float) alignment - # - # 'figure_align': 'htbp', -} - -# Grouping the document tree into LaTeX files. List of tuples -# (source start file, target name, title, -# author, documentclass [howto, manual, or own class]). +htmlhelp_basename = f'{project}doc' +doc_title = f'{project} Documentation' latex_documents = [ - (master_doc, 'diffxpy.tex', 'diffxpy Documentation', - 'Florian R. Hölzlwimmer, David S. Fischer', 'manual'), + (master_doc, f'{project}.tex', doc_title, author, 'manual'), ] - -# -- Options for manual page output ------------------------------------------ - -# One entry per manual page. List of tuples -# (source start file, name, description, authors, manual section). man_pages = [ - (master_doc, 'diffxpy', 'diffxpy Documentation', - [author], 1) + (master_doc, project, doc_title, [author], 1) ] - -# -- Options for Texinfo output ---------------------------------------------- - -# Grouping the document tree into Texinfo files. List of tuples -# (source start file, target name, title, author, -# dir menu entry, description, category) texinfo_documents = [ - (master_doc, 'diffxpy', 'diffxpy Documentation', - author, 'diffxpy', 'One line description of project.', - 'Miscellaneous'), + (master_doc, project, doc_title, author, project, 'One line description of project.', 'Miscellaneous'), ] -# -- Extension configuration ------------------------------------------------- - - -# -- generate_options override ------------------------------------------ -# TODO: why? - - -def process_generate_options(app: Sphinx): - genfiles = app.config.autosummary_generate - - if genfiles and not hasattr(genfiles, '__len__'): - env = app.builder.env - genfiles = [ - env.doc2path(x, base=None) - for x in env.found_docs - if Path(env.doc2path(x)).is_file() - ] - if not genfiles: - return - - from sphinx.ext.autosummary.generate import generate_autosummary_docs +# -- Images for plot functions ------------------------------------------------- - ext = app.config.source_suffix - genfiles = [ - genfile + (not genfile.endswith(tuple(ext)) and ext[0] or '') - for genfile in genfiles - ] - suffix = autosummary.get_rst_suffix(app) - if suffix is None: - return +def insert_function_images(app, what, name, obj, options, lines): + path = Path(__file__).parent / 'api' / f'{name}.png' + if what != 'function' or not path.is_file(): return + lines[0:0] = [f'.. image:: {path.name}', ' :width: 200', ' :align: right', ''] - generate_autosummary_docs( - genfiles, builder=app.builder, - warn=logger.warning, info=logger.info, - suffix=suffix, base_path=app.srcdir, - imported_members=True, app=app, - ) +# -- GitHub links -------------------------------------------------------------- -autosummary.process_generate_options = process_generate_options +def autolink(url_template, title_template='{}'): + from docutils import nodes -# -- GitHub URLs for class and method pages ------------------------------------------ - - -def get_obj_module(qualname): - """Get a module/class/attribute and its original module by qualname""" - modname = qualname - classname = None - attrname = None - while modname not in sys.modules: - attrname = classname - modname, classname = modname.rsplit('.', 1) - - # retrieve object and find original module name - if classname: - cls = getattr(sys.modules[modname], classname) - modname = cls.__module__ - obj = getattr(cls, attrname) if attrname else cls - else: - obj = None - - return obj, sys.modules[modname] - - -def get_linenos(obj): - """Get an object’s line numbers""" - try: - lines, start = inspect.getsourcelines(obj) - except TypeError: - return None, None - else: - return start, start + len(lines) - 1 - - -project_dir = Path(__file__).parent.parent # project/docs/conf.py/../.. → project/ -github_url1 = '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/{github_user}/{github_repo}/tree/{github_version}'.format_map(html_context) -github_url2 = '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/theislab/diffxpy/tree/master' - - -def modurl(qualname: str) -> str: - """Get the full GitHub URL for some object’s qualname.""" - obj, module = get_obj_module(qualname) - github_url = github_url1 - try: - path = Path(module.__file__).relative_to(project_dir) - except ValueError: - # trying to document something from another package - github_url = github_url2 - path = '/'.join(module.__file__.split('/')[-2:]) - start, end = get_linenos(obj) - fragment = f'#L{start}-L{end}' if start and end else '' - return f'{github_url}/{path}{fragment}' - - -def api_image(qualname: str) -> Optional[str]: - # I’d like to make this a contextfilter, but the jinja context doesn’t contain the path, - # so no chance to not hardcode “api/” here. - path = Path(__file__).parent / 'api' / f'{qualname}.png' - print(path, path.is_file()) - return f'.. image:: {path.name}\n :width: 200\n :align: right' if path.is_file() else '' - - -# html_context doesn’t apply to autosummary templates ☹ -# and there’s no way to insert filters into those templates -# so we have to modify the default filters -from jinja2.defaults import DEFAULT_FILTERS - -DEFAULT_FILTERS.update(modurl=modurl, api_image=api_image) - -# -- Prettier Param docs -------------------------------------------- - - -from typing import Dict, List, Tuple - -from docutils import nodes -from sphinx import addnodes -from sphinx.domains.python import PyTypedField, PyObject -from sphinx.environment import BuildEnvironment - - -class PrettyTypedField(PyTypedField): - list_type = nodes.definition_list - - def make_field( - self, - types: Dict[str, List[nodes.Node]], - domain: str, - items: Tuple[str, List[nodes.inline]], - env: BuildEnvironment = None - ) -> nodes.field: - def makerefs(rolename, name, node): - return self.make_xrefs(rolename, domain, name, node, env=env) - - def handle_item(fieldarg: str, content: List[nodes.inline]) -> nodes.definition_list_item: - head = nodes.term() - head += makerefs(self.rolename, fieldarg, addnodes.literal_strong) - fieldtype = types.pop(fieldarg, None) - if fieldtype is not None: - head += nodes.Text(' : ') - if len(fieldtype) == 1 and isinstance(fieldtype[0], nodes.Text): - typename = ''.join(n.astext() for n in fieldtype) - head += makerefs(self.typerolename, typename, addnodes.literal_emphasis) - else: - head += fieldtype - - body_content = nodes.paragraph('', '', *content) - body = nodes.definition('', body_content) - - return nodes.definition_list_item('', head, body) - - fieldname = nodes.field_name('', self.label) - if len(items) == 1 and self.can_collapse: - fieldarg, content = items[0] - bodynode = handle_item(fieldarg, content) - else: - bodynode = self.list_type() - for fieldarg, content in items: - bodynode += handle_item(fieldarg, content) - fieldbody = nodes.field_body('', bodynode) - return nodes.field('', fieldname, fieldbody) - - -# replace matching field types with ours -PyObject.doc_field_types = [ - PrettyTypedField( - ft.name, - names=ft.names, - typenames=ft.typenames, - label=ft.label, - rolename=ft.rolename, - typerolename=ft.typerolename, - can_collapse=ft.can_collapse, - ) if isinstance(ft, PyTypedField) else ft - for ft in PyObject.doc_field_types -] + def role(name, rawtext, text, lineno, inliner, options={}, content=[]): + url = url_template.format(text) + title = title_template.format(text) + node = nodes.reference(rawtext, title, refuri=url, **options) + return [node], [] + return role diff --git a/docs/requires.txt b/docs/requires.txt index 7d2d66e..228129e 100644 --- a/docs/requires.txt +++ b/docs/requires.txt @@ -2,7 +2,7 @@ numpy>=1.14.0 scipy pandas patsy>=0.5.0 -batchglm>=0.4.0 +batchglm>=0.6.0 xarray statsmodels anndata diff --git a/setup.py b/setup.py index 0a50c62..cb70fda 100644 --- a/setup.py +++ b/setup.py @@ -21,7 +21,7 @@ 'scipy', 'pandas', 'patsy>=0.5.0', - 'batchglm>=0.4.0', + 'batchglm>=0.6.0', 'xarray', 'statsmodels', ], @@ -29,10 +29,6 @@ 'optional': [ 'anndata', ], - # 'scanpy_deps': [ - # "scanpy", - # "anndata" - # ], 'plotting_deps': [ "seaborn", "matplotlib"