differential expression improvements (#452)

* add cutoff for low expression genes in topN selection

* remove debugging printfs

* change cli param name for diffexp cutoff

* change CLI param name

* second try at diffexp - using lfc sort with pval cutoff

* use lfc cutoff

* update comments to match code; cap p-value adjustment to max of 1

* lint

* explain diffexp in readme

* add link

* add diffexp-lfc-cutoff to test config

* update test to match revised diffexp spec

* fix latent bug in GET arg handling that was breaking tests

* lint

* comment cleanup

* fix variance overestimation so it is symmetric

* lint

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.

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_
 

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

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'])

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]

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):
         """

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
     }

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)