diff --git a/README.md b/README.md index 518fbaef..6ec6e421 100644 --- a/README.md +++ b/README.md @@ -151,6 +151,10 @@ Currently this is not supported directly, but you should be able to do this manu This may happen, especially as we work out bugs in our installation process! Please create a new [Github issue](https://github.com/chanzuckerberg/cellxgene/issues), explain what you did, and include all the error messages you saw. It'd also be super helpful if you call `pip freeze` and include the full output alongside your issue. +> How are you computing and sorting differential expression results? + +Currently we use a [Welch's *t*-test](https://en.wikipedia.org/wiki/Welch%27s_t-test) implementation including the same variance overestimation correction as used in `scanpy`. We sort the `tscore` to identify the top N genes, and then filter to remove any that fall below a cutoff log fold change value, which can help remove spurious test results. The default threshold is `0.01` and can be changed using the option `--diffexp-lfc-cutoff`. We can explore adding support for other test types in the future. + > I'm following the developer instructions and get an error about "missing files and directories” when trying to build the client This is likely because you do not have node and npm installed, we recommend using [nvm](https://github.com/creationix/nvm) if you're new to using these tools. diff --git a/docs/REST_API.md b/docs/REST_API.md index 9ebc71b6..3313698b 100644 --- a/docs/REST_API.md +++ b/docs/REST_API.md @@ -568,14 +568,14 @@ If differential expression is not supported by the server, must return an HTTP 5 **Response body:** -- For 200 Success, differential expression statistics returned as array of arrays sorted by varindex, where each contains the following values: +- For 200 Success, differential expression statistics returned as array of arrays, where each contains the following values: - **varIndex**: variable index for the computed results - **logfoldchange**: log fold-change of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group, - **pVal**: unadjusted p-value, - **pValAdj**: adjusted p-value - Statistics are encoded as an array of arrays, with fields ordered as: + Values ordered as: _varIndex_, _logfoldchange_, _pVal_, _pValAdj_ diff --git a/server/app/driver/driver.py b/server/app/driver/driver.py index b4693336..fd5c352d 100644 --- a/server/app/driver/driver.py +++ b/server/app/driver/driver.py @@ -17,6 +17,7 @@ class CXGDriver(metaclass=ABCMeta): self.layout_method = args["layout"] self.diffexp_method = args["diffexp"] self.max_category_items = args["max_category_items"] + self.diffexp_lfc_cutoff = args["diffexp_lfc_cutoff"] self.cluster = None @property diff --git a/server/app/rest_api/rest.py b/server/app/rest_api/rest.py index fecd7dd6..7e613c5b 100644 --- a/server/app/rest_api/rest.py +++ b/server/app/rest_api/rest.py @@ -346,7 +346,7 @@ class DataObsAPI(Resource): def get(self): accept_type = request.args.get("accept-type", None) # request.args is immutable - args = dict(request.args) + args = request.args.copy() args.pop("accept-type", None) try: filter_ = parse_filter(ImmutableMultiDict(args), current_app.data.schema['annotations']) @@ -450,7 +450,7 @@ class DataVarAPI(Resource): def get(self): accept_type = request.args.get("accept-type", None) # request.args is immutable - args = dict(request.args) + args = request.args.copy() args.pop("accept-type", None) try: filter_ = parse_filter(ImmutableMultiDict(args), current_app.data.schema['annotations']) diff --git a/server/app/scanpy_engine/diffexp.py b/server/app/scanpy_engine/diffexp.py index d9c0e2a4..8f229b7d 100644 --- a/server/app/scanpy_engine/diffexp.py +++ b/server/app/scanpy_engine/diffexp.py @@ -25,25 +25,42 @@ def _mean_var_n(X): return mean, v, n -def diffexp_ttest(adata, maskA, maskB, top_n=8): +def diffexp_ttest(adata, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01): """ - Return differential expression statistics for top N variables, sorted by - t statistic. Implemented as a unequal variance t-test. + Return differential expression statistics for top N variables. + + Algorithm: + - compute log fold change (log2(meanA/meanB)) + - compute Welch's t-test statistic and pvalue (w/ Bonferroni correction) + - return top N abs(logfoldchange) where lfc > diffexp_lfc_cutoff + + If there are not N which meet criteria, augment by removing the logfoldchange + threshold requirement. + + Notes on alogrithm: + - Welch's ttest provides basic statistics test. + https://en.wikipedia.org/wiki/Welch%27s_t-test + - p-values adjusted with Bonferroni correction. + https://en.wikipedia.org/wiki/Bonferroni_correction :param adata: anndata dataframe :param maskA: observation selection mask for set 1 :param maskB: observation selection mask for set 2 :param top_n: number of variables to return stats for + :param diffexp_lfc_cutoff: minimum :return: for top N genes, [ varindex, logfoldchange, pval, pval_adj ] """ - # mean, variance, N + if top_n > adata.n_obs: + top_n = adata.n_obs + + # mean, variance, N - calculate for both selections meanA, vA, nA = _mean_var_n(adata._X[maskA]) meanB, vB, nB = _mean_var_n(adata._X[maskB]) # variance / N - vnA = vA / nA - vnB = vB / nB + vnA = vA / min(nA, nB) # overestimate variance, would normally be nA + vnB = vB / min(nA, nB) # overestimate variance, would normally be nB sum_vn = vnA + vnB # degrees of freedom for Welch's t-test @@ -59,18 +76,31 @@ def diffexp_ttest(adata, maskA, maskB, top_n=8): # p-value pvals = stats.t.sf(np.abs(tscores), dof) * 2 pvals_adj = pvals * adata._X.shape[1] + pvals_adj[pvals_adj > 1] = 1 # cap adjusted p-value at 1 # logfoldchanges: log2(meanA / meanB) logfoldchanges = np.log2(np.abs((meanA + 1e-9) / (meanB + 1e-9))) - # top n sort + # find all with lfc > cutoff + lfc_above_cutoff_idx = np.nonzero(np.abs(logfoldchanges) > diffexp_lfc_cutoff)[0] stats_to_sort = np.abs(tscores) - partition = np.argpartition(stats_to_sort, -top_n)[-top_n:] - rel_sort_order = np.argsort(stats_to_sort[partition])[::-1] - vars_indices = np.arange(adata.n_vars, dtype=int) - sort_order = vars_indices[partition][rel_sort_order] - # top n slice + # derive sort order + if lfc_above_cutoff_idx.shape[0] > top_n: + # partition top N + rel_t_partition = np.argpartition(stats_to_sort[lfc_above_cutoff_idx], -top_n)[-top_n:] + t_partition = lfc_above_cutoff_idx[rel_t_partition] + # sort the top N partition + rel_sort_order = np.argsort(stats_to_sort[t_partition])[::-1] + sort_order = t_partition[rel_sort_order] + else: + # partition and sort top N, ignoring lfc cutoff + partition = np.argpartition(stats_to_sort, -top_n)[-top_n:] + rel_sort_order = np.argsort(stats_to_sort[partition])[::-1] + indices = np.indices(stats_to_sort.shape)[0] + sort_order = indices[partition][rel_sort_order] + + # top n slice based upon sort order logfoldchanges_top_n = logfoldchanges[sort_order] pvals_top_n = pvals[sort_order] pvals_adj_top_n = pvals_adj[sort_order] diff --git a/server/app/scanpy_engine/scanpy_engine.py b/server/app/scanpy_engine/scanpy_engine.py index 79dbe328..8b1845a2 100644 --- a/server/app/scanpy_engine/scanpy_engine.py +++ b/server/app/scanpy_engine/scanpy_engine.py @@ -325,8 +325,8 @@ class ScanpyEngine(CXGDriver): raise FilterError(f"Error parsing filter: {e}") from e if top_n is None: top_n = DEFAULT_TOP_N - result = diffexp_ttest(self.data, obs_mask_A, obs_mask_B, top_n) - return sorted(result, key=lambda r: r[0]) + result = diffexp_ttest(self.data, obs_mask_A, obs_mask_B, top_n, self.diffexp_lfc_cutoff) + return result def layout(self, filter, interactive_limit=None): """ diff --git a/server/cli/launch.py b/server/cli/launch.py index 975f4132..ab5396c1 100644 --- a/server/cli/launch.py +++ b/server/cli/launch.py @@ -27,8 +27,10 @@ from server.app.util.errors import ScanpyFileError help="Bind to all interfaces (this makes the server accessible beyond this computer).") @click.option("--max-category-items", default=100, metavar="", show_default=True, help="Limits the number of categorical annotation items displayed.") +@click.option("--diffexp-lfc-cutoff", default=0.01, show_default=True, + help="Relative expression cutoff used when selecting top N differentially expressed genes") def launch(data, layout, diffexp, title, verbose, debug, obs_names, var_names, - open_browser, port, listen_all, max_category_items): + open_browser, port, listen_all, max_category_items, diffexp_lfc_cutoff): """Launch the cellxgene data viewer. This web app lets you explore single-cell expression data. Data must be in a format that cellxgene expects, read the @@ -92,6 +94,7 @@ def launch(data, layout, diffexp, title, verbose, debug, obs_names, var_names, "layout": layout, "diffexp": diffexp, "max_category_items": max_category_items, + "diffexp_lfc_cutoff": diffexp_lfc_cutoff, "obs_names": obs_names, "var_names": var_names } diff --git a/server/test/test_scanpy_engine.py b/server/test/test_scanpy_engine.py index 7ab0b971..0b8ccdcd 100644 --- a/server/test/test_scanpy_engine.py +++ b/server/test/test_scanpy_engine.py @@ -14,7 +14,7 @@ from server.app.scanpy_engine.scanpy_engine import ScanpyEngine class UtilTest(unittest.TestCase): def setUp(self): args = {'layout': 'umap', 'diffexp': 'ttest', 'max_category_items': 100, - 'obs_names': None, 'var_names': None} + 'obs_names': None, 'var_names': None, 'diffexp_lfc_cutoff': 0.01} self.data = ScanpyEngine("example-dataset/pbmc3k.h5ad", args) self.data._create_schema() @@ -201,8 +201,6 @@ class UtilTest(unittest.TestCase): } result = self.data.diffexp_topN(f1["filter"], f2["filter"]) self.assertEqual(len(result), 10) - var_idx = [i[0] for i in result] - self.assertEqual(var_idx, sorted(var_idx)) result = self.data.diffexp_topN(f1["filter"], f2["filter"], 20) self.assertEqual(len(result), 20)