Add support for anndata backed mode (#943)

* initial cut at backed mode

* make flask multithreading conditional on debug flag

* update X access to support backed mode

* lint

* improve help message for backed mode

* fix tests

* add MatrixProxy to normalize supported matrix types

* add FAQ entry for --backed

* remove use of matrix.T

* clean up

* add ability to disable diffexp from CLI; add hueristic to detect likely slow diffexp calculation, and warn user

* fix tests

* do not print diffexp speed warning if diffexp is disabled

* tweak wording of diffexp speed messages

* add FAQ entry on --disable-diffexp

* revise heuristic for warning about slow diffexp

* use quick tooltip delay on diffexp button

client/src/components/menubar/index.js

@@ -16,7 +16,10 @@ import CellSetButton from "./cellSetButtons";
 import InformationMenu from "./infoMenu";
 import UndoRedoReset from "./undoRedoReset";
 import Clip from "./clip";
-import { tooltipHoverOpenDelay } from "../../globals";
+import {
+  tooltipHoverOpenDelay,
+  tooltipHoverOpenDelayQuick
+} from "../../globals";
 
 @connect(state => ({
   universe: state.universe,
@@ -39,7 +42,9 @@ import { tooltipHoverOpenDelay } from "../../globals";
   libraryVersions: state.config?.library_versions, // eslint-disable-line camelcase
   undoDisabled: state["@@undoable/past"].length === 0,
   redoDisabled: state["@@undoable/future"].length === 0,
-  aboutLink: state.config?.links?.["about-dataset"]
+  aboutLink: state.config?.links?.["about-dataset"],
+  disableDiffexp: state.config?.parameters?.["disable-diffexp"] ?? false,
+  diffexpMayBeSlow: state.config?.parameters?.["diffexp-may-be-slow"] ?? false
 }))
 class MenuBar extends React.Component {
   static isValidDigitKeyEvent(e) {
@@ -255,6 +260,65 @@ class MenuBar extends React.Component {
     });
   };
 
+  renderDiffExp() {
+    /* diffexp-related buttons may be disabled */
+    const { disableDiffexp, differential, diffexpMayBeSlow } = this.props;
+    if (disableDiffexp) return null;
+
+    const haveBothCellSets =
+      !!differential.celllist1 && !!differential.celllist2;
+
+    const tipMessage =
+      "See top 10 differentially expressed genes" +
+      (diffexpMayBeSlow
+        ? " (CAUTION: large dataset - may take longer or fail)"
+        : "");
+
+    return (
+      <div className="bp3-button-group" style={{ marginRight: 10 }}>
+        <CellSetButton {...this.props} eitherCellSetOneOrTwo={1} />
+        <CellSetButton {...this.props} eitherCellSetOneOrTwo={2} />
+        {!differential.diffExp ? (
+          <Tooltip
+            content={tipMessage}
+            position="bottom"
+            hoverOpenDelay={tooltipHoverOpenDelayQuick}
+          >
+            <AnchorButton
+              disabled={!haveBothCellSets}
+              intent="primary"
+              data-testid="diffexp-button"
+              loading={differential.loading}
+              icon="left-join"
+              fill
+              type="button"
+              onClick={this.computeDiffExp}
+            >
+              Compute Differential Expression
+            </AnchorButton>
+          </Tooltip>
+        ) : null}
+
+        {differential.diffExp ? (
+          <Tooltip
+            content="Remove differentially expressed gene list and clear cell selections"
+            position="bottom"
+            hoverOpenDelay={tooltipHoverOpenDelayQuick}
+          >
+            <Button
+              type="button"
+              fill
+              intent="warning"
+              onClick={this.clearDifferentialExpression}
+            >
+              Clear Differential Expression
+            </Button>
+          </Tooltip>
+        ) : null}
+      </div>
+    );
+  }
+
   render() {
     const {
       dispatch,
@@ -273,9 +337,6 @@ class MenuBar extends React.Component {
     } = this.props;
     const { pendingClipPercentiles } = this.state;
 
-    const haveBothCellSets =
-      !!differential.celllist1 && !!differential.celllist2;
-
     // constants used to create selection tool button
     let selectionTooltip;
     let selectionButtonClass;
@@ -295,47 +356,7 @@ class MenuBar extends React.Component {
           top: 8
         }}
       >
-        <div className="bp3-button-group" style={{ marginRight: 10 }}>
+        {this.renderDiffExp()}
-          <CellSetButton {...this.props} eitherCellSetOneOrTwo={1} />
-          <CellSetButton {...this.props} eitherCellSetOneOrTwo={2} />
-          {!differential.diffExp ? (
-            <Tooltip
-              content="Add two cells selections, see the top 15 differentially expressed genes between them"
-              position="bottom"
-              hoverOpenDelay={tooltipHoverOpenDelay}
-            >
-              <AnchorButton
-                disabled={!haveBothCellSets}
-                intent="primary"
-                data-testid="diffexp-button"
-                loading={differential.loading}
-                icon="left-join"
-                fill
-                type="button"
-                onClick={this.computeDiffExp}
-              >
-                Compute Differential Expression
-              </AnchorButton>
-            </Tooltip>
-          ) : null}
-
-          {differential.diffExp ? (
-            <Tooltip
-              content="Remove differentially expressed gene list and clear cell selections"
-              position="bottom"
-              hoverOpenDelay={tooltipHoverOpenDelay}
-            >
-              <Button
-                type="button"
-                fill
-                intent="warning"
-                onClick={this.clearDifferentialExpression}
-              >
-                Clear Differential Expression
-              </Button>
-            </Tooltip>
-          ) : null}
-        </div>
         <Tooltip
           content="Show only metadata and cells which are currently selected"
           position="bottom"

client/src/globals.js

@@ -14,7 +14,9 @@ export const configDefaults = {
   features: {},
   displayNames: {},
   parameters: {
-    "max-category-items": 1000
+    "max-category-items": 1000,
+    "disable-diffexp": false,
+    "diffexp-may-be-slow": false
   },
   links: {}
 };

docs/faq.md

@@ -12,7 +12,7 @@ Currently, you can go straight into `cellxgene launch` with your own analyzed da
 
 If your data is in a different format, and/or you still need to perform dimensionality reduction and clustering, `cellxgene` can do that for you with the `prepare` command. `cellxgene prepare` runs `scanpy` under the hood and can read in any format that is currently supported by `scanpy` (including mtx, loom, and more listed [here](https://scanpy.readthedocs.io/en/latest/api/index.html#reading)).
 
-To add `cellxgene prepare` to your cellxgene installation run `pip install cellxgene[prepare]`. 
+To add `cellxgene prepare` to your cellxgene installation run `pip install cellxgene[prepare]`.
 
 The output of `cellxgene prepare` is a h5ad file with your computed clusters and tsne/umap projections that can be used in `cellxgene launch`.
 
@@ -110,3 +110,21 @@ For example:
 pip install s3fs
 cellxgene launch s3://mybucket.s3-us-west-2.amazonaws.com/mydata.h5ad
 ```
+
+#### What does the command line option `--backed` do?
+
+The `--backed` option instructs `cellxgene launch` to read the H5AD file in "backed" mode (for more information, see the
+[anndata.read_h5ad() documentation](https://anndata.readthedocs.io/en/latest/anndata.read_h5ad.html#anndata.read_h5ad)).
+
+By default, cellxgene will read the entire H5AD will be into memory at startup, improving application speed and performance.
+Very large datasets may not fit in memory. The "--backed" mode will read the file incrementally, reducing memory
+use, and for large files, improving startup speed. _However_, this option will also significantly slow
+down access to gene expression histograms, and may render differential expression calculations too slow
+to use (see `--disable-diffexp` for an option to disable this feature).
+
+#### What does the command line option `--disable-diffexp` do?
+
+The `--disable-diffexp` option will disable and hide the `Compute Differential Expression` feature.
+For large datasets, or datasets loaded with the `--backed` option, computing differential expression may
+be extremely slow or use excessive reources on the host computer (eg, memory thrasing).
+Disabling the feature will ensure that the end-user does not accidentally initiate this computation.

server/app/driver/driver.py

@@ -29,7 +29,9 @@ class CXGDriver(metaclass=ABCMeta):
         return {
             "layout": None,
             "max_category_items": None,
-            "diffexp_lfc_cutoff": None
+            "diffexp_lfc_cutoff": None,
+            "disable_diffexp": False,
+            "diffexp_may_be_slow": False
         }
 
     @property

server/app/rest_api/rest.py

@@ -69,7 +69,9 @@ class ConfigAPI(Resource):
                     "about-dataset": current_app.config["ABOUT_DATASET"]
                 },
                 "parameters": {
-                    "max-category-items": current_app.data.config["max_category_items"]
+                    "max-category-items": current_app.data.config["max_category_items"],
+                    "disable-diffexp": current_app.data.config["disable_diffexp"],
+                    "diffexp-may-be-slow": current_app.data.config["diffexp_may_be_slow"]
                 },
                 "library_versions": {
                     "cellxgene": cellxgene_version,

server/app/scanpy_engine/diffexp.py

@@ -62,13 +62,12 @@ def diffexp_ttest(adata, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
     :param diffexp_lfc_cutoff: minimum
     :return:  for top N genes, [ varindex, logfoldchange, pval, pval_adj ]
     """
-
     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])
+    meanA, vA, nA = _mean_var_n(adata.X[maskA, :])
-    meanB, vB, nB = _mean_var_n(adata._X[maskB])
+    meanB, vB, nB = _mean_var_n(adata.X[maskB, :])
 
     # variance / N
     vnA = vA / min(nA, nB)  # overestimate variance, would normally be nA
@@ -87,7 +86,7 @@ def diffexp_ttest(adata, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
 
     # p-value
     pvals = stats.t.sf(np.abs(tscores), dof) * 2
-    pvals_adj = pvals * adata._X.shape[1]
+    pvals_adj = pvals * adata.X.shape[1]
     pvals_adj[pvals_adj > 1] = 1  # cap adjusted p-value at 1
 
     # logfoldchanges: log2(meanA / meanB)

server/app/scanpy_engine/matrix_proxy.py

@@ -0,0 +1,42 @@
+
+from server.app.util.matrix_proxy import MatrixProxyView, ArrayProxyView
+
+"""
+AnnData/h5py are inconsistent in the API supported by various types of
+X matrices.  Sometimes you get a fully ndarray, sometims a Scipy sparse
+matrix, sometimes h5py proxies with a subset of our needed interfaces.
+
+This glue code paves over all of that, providing the core set of methods
+that the cellxgene ScanPy driver assumes are in existance.  Put another
+way, all of the non-portable assumptions are here.
+"""
+
+
+class ArrayProxyView_anndata_h5py(ArrayProxyView):
+    """
+    override to handle sparse getitem semantics, which differ
+    from numpy.
+    """
+    def toarray(self):
+        """ sadly, sparse indexing doesn't drop dimensions like numpy! """
+        arr = self.m[self._index[0], self._index[1]]
+        if self._vdim == 0:
+            arr = arr.transpose()
+        return arr.toarray()[0]
+
+
+class MatrixProxy_anndata_h5py(MatrixProxyView):
+    """
+    AnnData sparse array stored in H5AD, or proxies for backed data.
+    None of these handle indexing very well, so we plop a proxy on top.
+    """
+    @classmethod
+    def __supports__(cls):
+        return ("anndata.h5py.h5sparse.SparseDataset",
+                "anndata.h5py.h5sparse.backed_csc_matrix",
+                "anndata.h5py.h5sparse.backed_csr_matrix",
+                "h5py._hl.dataset.Dataset")
+
+    @classmethod
+    def create_array(cls, *args, **kwargs):
+        return ArrayProxyView_anndata_h5py(*args, **kwargs)

server/app/scanpy_engine/scanpy_engine.py

@@ -21,6 +21,14 @@ from server.app.util.utils import jsonify_scanpy, requires_data
 from server.app.scanpy_engine.diffexp import diffexp_ttest
 from server.app.util.fbs.matrix import encode_matrix_fbs, decode_matrix_fbs
 from server.app.scanpy_engine.labels import read_labels, write_labels
+import server.app.scanpy_engine.matrix_proxy  # noqa: F401
+from server.app.util.matrix_proxy import MatrixProxy
+
+
+def has_method(o, name):
+    """ return True if `o` has callable method `name` """
+    op = getattr(o, name, None)
+    return op is not None and callable(op)
 
 
 class ScanpyEngine(CXGDriver):
@@ -45,6 +53,9 @@ class ScanpyEngine(CXGDriver):
             "var_names": None,
             "diffexp_lfc_cutoff": 0.01,
             "label_file": None,
+            "backed": False,
+            "disable_diffexp": False,
+            "diffexp_may_be_slow": False
         }
 
     @staticmethod
@@ -208,7 +219,8 @@ class ScanpyEngine(CXGDriver):
             with data_locator.local_handle() as lh:
                 # as of AnnData 0.6.19, backed mode performs initial load fast, but at the
                 # cost of significantly slower access to X data.
-                self.data = anndata.read_h5ad(lh)
+                backed = 'r' if self.config['backed'] else None
+                self.data = anndata.read_h5ad(lh, backed=backed)
 
         except ValueError:
             raise ScanpyFileError(
@@ -251,6 +263,11 @@ class ScanpyEngine(CXGDriver):
         self._validate_label_data()
         self._create_schema()
 
+        # heuristic
+        n_values = self.data.shape[0] * self.data.shape[1]
+        if (n_values > 1e8 and self.config['backed'] is True) or (n_values > 5e8):
+            self.config.update({"diffexp_may_be_slow": True})
+
     @requires_data
     def _default_and_validate_layouts(self):
         """ function:
@@ -467,22 +484,6 @@ class ScanpyEngine(CXGDriver):
 
         return jsonify_scanpy({"status": "OK"})
 
-    @staticmethod
-    def slice_columns(X, var_mask):
-        """
-        Slice columns from the matrix X, as specified by the mask
-        Semantically equivalent to X[:, var_mask], but handles sparse
-        matrices in a more performant manner.
-        """
-        if var_mask is None:    # noop
-            return X
-        if sparse.issparse(X):  # use tuned getcol/hstack for performance
-            indices = np.nonzero(var_mask)[0]
-            cols = [X.getcol(i) for i in indices]
-            return sparse.hstack(cols, format="csc")
-        else:   # else, just use standard slicing, which is fine for dense arrays
-            return X[:, var_mask]
-
     @requires_data
     def data_frame_to_fbs_matrix(self, filter, axis):
         """
@@ -505,7 +506,8 @@ class ScanpyEngine(CXGDriver):
             raise FilterError("filtering on obs unsupported")
 
         # Currently only handles VAR dimension
-        X = self.slice_columns(self.data._X, var_selector)
+        X = MatrixProxy.create(self.data.X if var_selector is None
+                               else self.data.X[:, var_selector])
         return encode_matrix_fbs(X, col_idx=np.nonzero(var_selector)[0], row_idx=None)
 
     @requires_data

server/app/util/fbs/matrix.py

@@ -12,6 +12,7 @@ import server.app.util.fbs.NetEncoding.Uint32Array as Uint32Array
 import server.app.util.fbs.NetEncoding.Float32Array as Float32Array
 import server.app.util.fbs.NetEncoding.Float64Array as Float64Array
 import server.app.util.fbs.NetEncoding.JSONEncodedArray as JSONEncodedArray
+from server.app.util.matrix_proxy import MatrixProxy
 
 
 # Placeholder until recent enhancements to flatbuffers Python
@@ -24,7 +25,7 @@ def CreateNumpyVector(builder, x):
     """CreateNumpyVector writes a numpy array into the buffer."""
 
     if not isinstance(x, np.ndarray):
-        raise TypeError("non-numpy-ndarray passed to CreateNumpyVector")
+        raise TypeError(f"non-numpy-ndarray passed to CreateNumpyVector ({type(x)}")
 
     if x.dtype.kind not in ['b', 'i', 'u', 'f']:
         raise TypeError("numpy-ndarray holds elements of unsupported datatype")
@@ -91,7 +92,7 @@ def serialize_typed_array(builder, source_array, encoding_info):
         as_json = arr.to_json(orient='records')
         arr = np.array(bytearray(as_json, 'utf-8'))
     else:
-        if sparse.issparse(arr):
+        if MatrixProxy.ismatrixproxy(arr) or sparse.issparse(arr):
             arr = arr.toarray()
         elif isinstance(arr, pd.Series):
             arr = arr.get_values()
@@ -99,8 +100,11 @@ def serialize_typed_array(builder, source_array, encoding_info):
             arr = arr.astype(as_type)
 
     # serialize the ndarray into a vector
-    if arr.ndim == 2 and arr.shape[0] == 1:
+    if arr.ndim == 2:
-        arr = arr[0]
+        if arr.shape[0] == 1:
+            arr = arr[0]
+        elif arr.shape[1] == 1:
+            arr = arr.T[0]
     vec = CreateNumpyVector(builder, arr)
 
     # serialize the typed array table
@@ -185,16 +189,11 @@ def encode_matrix_fbs(matrix, row_idx=None, col_idx=None):
     # estimate size needed, so we don't unnecessarily realloc.
     builder = flatbuffers.Builder(guess_at_mem_needed(matrix))
 
-    if isinstance(matrix, pd.DataFrame):
-        matrix_columns = reversed(tuple(matrix[name] for name in matrix))
-    else:
-        matrix_columns = reversed(tuple(c for c in matrix.T))
-
     columns = []
-    # for idx in reversed(np.arange(n_cols)):
+    for cidx in range(n_cols - 1, -1, -1):
-    for c in matrix_columns:
         # serialize the typed array
-        typed_arr = serialize_typed_array(builder, c, column_encoding)
+        col = matrix.iloc[:, cidx] if isinstance(matrix, pd.DataFrame) else matrix[:, cidx]
+        typed_arr = serialize_typed_array(builder, col, column_encoding)
 
         # serialize the Column union
         columns.append(serialize_column(builder, typed_arr))

server/app/util/matrix_proxy.py

@@ -0,0 +1,425 @@
+import abc
+from itertools import zip_longest
+from copy import copy
+import numpy as np
+
+"""
+cellxgene deals with a variety of matrix data types, many of which do
+not support a consistent API.  This framework allows proxies to be created
+to pave over some of this.  Most significantly, AnnData.X does not guarantee
+that much of its API (eg. .X.T) will work.
+"""
+
+INT_TYPES = (int, np.integer)
+
+
+class _ArrayProxyBase(abc.ABC):
+    """
+    Private base class for array or matrix proxy.  This summarizes
+    the interface used by the rest of cellxgene.
+    """
+    @property
+    @abc.abstractmethod
+    def dtype(self):
+        raise NotImplementedError()
+
+    @property
+    @abc.abstractmethod
+    def ndim(self):
+        raise NotImplementedError()
+
+    @property
+    @abc.abstractmethod
+    def shape(self):
+        raise NotImplementedError()
+
+    @property
+    @abc.abstractmethod
+    def T(self):
+        raise NotImplementedError()
+
+    @abc.abstractmethod
+    def __iter__(self):
+        raise NotImplementedError()
+
+    @abc.abstractmethod
+    def __getitem__(self, args):
+        raise NotImplementedError()
+
+    @abc.abstractmethod
+    def toarray():
+        raise NotImplementedError()
+
+
+class MatrixProxy(_ArrayProxyBase):
+    """
+    Abstract class - all interfaces we need, plus a factory method
+    to create a proxy based upon actual matrix type.
+
+    This class primarily provides the factory method and related support.
+    All other functionality is delegated to subclasses.
+    """
+
+    """
+    Registry of types to proxy class, where values are:
+    * None: unsupported
+    * True: self-supported
+    * string: proxy class
+    Sub-classes automatically register.
+    """
+    base_proxy_registry = {
+        'pandas.core.frame.DataFrame': True,
+        'numpy.ndarray': True,
+        'scipy.sparse.csc.csc_matrix': True,
+        'scipy.sparse.csr.csr_matrix': True,
+    }
+    proxy_registry = None
+    last_cache_token = None
+
+    @staticmethod
+    def _register_subclasses(subclasses, registry):
+        for c in subclasses:
+            names = c.__supports__()
+            for name in names:
+                registry[name] = c
+            MatrixProxy._register_subclasses(c.__subclasses__(), registry)
+
+    @classmethod
+    def build_proxy_registry(cls):
+        if cls.proxy_registry and abc.get_cache_token() == cls.last_cache_token:
+            return
+
+        cls.last_cache_token = abc.get_cache_token()
+        registry = copy(cls.base_proxy_registry)
+        MatrixProxy._register_subclasses(cls.__subclasses__(), registry)
+        cls.proxy_registry = registry
+
+    @classmethod
+    def create(cls, matrix):
+        """
+        Factory - call with a matrix and it will create a proxy if needed.
+        If the type already supports the necessary API, it is just returned
+        directly.
+        """
+        cls.build_proxy_registry()
+        t = type(matrix)
+        fqtn = t.__module__ + '.' + t.__name__
+        proxy_cls = cls.proxy_registry.get(fqtn, None)
+        if proxy_cls is None:
+            raise Exception(f"Matrix format `{fqtn}` is unsupported by proxy.")
+        if proxy_cls is True:
+            return matrix
+        return proxy_cls(matrix)
+
+    def __init__(self, m):
+        self.m = m
+
+    @classmethod
+    @abc.abstractmethod
+    def __supports__(cls):
+        raise NotImplementedError()
+
+    @classmethod
+    def ismatrixproxy(cls, m):
+        return isinstance(m, _ArrayProxyBase)
+
+
+class MatrixProxyView(MatrixProxy):
+    """
+    2D matrix view to a 2D matrix
+    """
+    def __init__(self, arg1, shape=None, index=(),
+                 transposed=False, copy=False):
+        if not copy:
+            m = arg1
+            super().__init__(m)
+
+            if shape is None:
+                shape = m.shape
+            assert(len(shape) == 2)
+
+            index = tuple(
+                map(lambda s_i:
+                    slice(0, s_i[0], 1) if s_i[1] is None else s_i[1],
+                    zip_longest(shape, index))
+            )
+
+            self._shape = shape
+            self._index = index
+            self.transposed = transposed
+
+        else:  # copy mode
+            super().__init__(arg1.m)
+            self._shape = arg1._shape
+            self._index = arg1._index
+            self.transposed = arg1.transposed
+
+    @classmethod
+    def create_array(cls, *args, **kwargs):
+        """ override if you use a different 1D array proxy """
+        return ArrayProxyView(*args, **kwargs)
+
+    def copy(self):
+        """ override if you need additional behaviors """
+        return self.__class__(self, copy=True)
+
+    @classmethod
+    def __supports__(cls):
+        return ()
+
+    def _swap(self, x):
+        return (x[1], x[0]) if self.transposed else x
+
+    @property
+    def dtype(self):
+        return self.m.dtype
+
+    @property
+    def ndim(self):
+        return len(self._shape)
+
+    @property
+    def shape(self):
+        return self._swap(self._shape)
+
+    @property
+    def T(self):
+        m = self.copy()
+        m.transposed = not m.transposed
+        return m
+
+    def __iter__(self):
+        M = self._swap(self._shape)[0]
+        s = self._swap((self._index))[0]
+        start, stop, step = s.indices(M)
+        for d in range(start, stop, step):
+            yield self[d]
+
+    def __getitem__(self, args):
+        """
+        decompose into the indexing patterns we use, throw for the rest.
+        Subclassses implement specialized access.
+        """
+        row, col = self._swap(_unpack_index(args, self.shape))
+        M, N = self._shape
+        allM, allN = self.m.shape
+
+        if isinstance(row, INT_TYPES):
+            row += self._index[0].start
+        elif isinstance(row, slice):
+            row = _slice_slice(self._index[0], allM, row, M)
+        if isinstance(col, INT_TYPES):
+            col += self._index[1].start
+        elif isinstance(col, slice):
+            col = _slice_slice(self._index[1], allN, col, N)
+
+        if isinstance(row, INT_TYPES):
+            if isinstance(col, INT_TYPES):
+                return self._getitem_intXint(row, col)
+            elif isinstance(col, slice):
+                return self._getitem_intXslice(row, col)
+        elif isinstance(row, slice):
+            if isinstance(col, INT_TYPES):
+                return self._getitem_sliceXint(row, col)
+            elif isinstance(col, slice):
+                return self._getitem_sliceXslice(row, col)
+
+        raise IndexError("unsupported column index types")
+
+    """
+    These getitem signatures are separate so that they may be
+    overridden by subclasses as necessary.  We don't do much
+    with them by default other than the obvious sub-slicing.
+
+    NOTE: these follow the numpy rules for dimensionality reduction
+    when an integer index is specified.
+    """
+    def _getitem_intXint(self, row, col):
+        return self.m[row, col]
+
+    def _getitem_intXslice(self, row, col):
+        shape = (_slice_length(col, self.m.shape[1]), )
+        return self.__class__.create_array(self.m, shape=shape, index=(row, col))
+
+    def _getitem_sliceXint(self, row, col):
+        shape = (_slice_length(row, self.m.shape[0]), )
+        return self.__class__.create_array(self.m, shape=shape, index=(row, col))
+
+    def _getitem_sliceXslice(self, row, col):
+        shape = (_slice_length(row, self.m.shape[0]),
+                 _slice_length(col, self.m.shape[1]))
+        return self.__class__(self.m, shape=shape, index=(row, col), transposed=self.transposed)
+
+    def toarray(self):
+        arr = self.m[self._index]
+        if self.transposed:
+            arr = arr.transpose()
+        return arr
+
+
+class ArrayProxyView(_ArrayProxyBase):
+    """
+    1D array view to a 2D matrix
+    """
+    def __init__(self, arg1, shape=None, index=None, copy=False):
+        super().__init__()
+        if not copy:
+            m = arg1
+
+            # one index MUST be an integer and the other MUST be a slice
+            assert(len(index) == 2)
+            assert(all(isinstance(idx, INT_TYPES + (slice, )) for idx in index))
+            assert(isinstance(index[0], INT_TYPES) != isinstance(index[1], INT_TYPES))
+
+            if shape is None:
+                if isinstance(index[0], INT_TYPES):
+                    shape = (m.shape[0], )
+                else:
+                    shape = (m.shape[1], )
+            assert(len(shape) == 1)
+
+            self._shape = shape
+            self.m = m
+            self._index = index
+            self._vdim = 1 if isinstance(index[0], INT_TYPES) else 0
+
+        else:
+            self.m = arg1.m
+            self._shape = arg1._shape
+            self._index = arg1._index
+            self.fixed = arg1.fixed
+
+    def copy(self):
+        return self.__class__(self, copy=True)
+
+    @property
+    def dtype(self):
+        return self.m.dtype
+
+    @property
+    def ndim(self):
+        return len(self._shape)
+
+    @property
+    def shape(self):
+        return self._shape
+
+    @property
+    def T(self):
+        return self.copy()
+
+    def __iter__(self):
+        _vdim = self._vdim
+        M = self._shape[0]
+        for d in range(*self._index[_vdim].indices(M)):
+            yield self[d]
+
+    def __getitem__(self, args):
+        _vdim = self._vdim
+        index = _unpack_index(args, self.shape)[0]
+        M = self._shape[0]
+        allM = self.m.shape[_vdim]
+
+        if isinstance(index, INT_TYPES):
+            index += self._index[_vdim].start
+        elif isinstance(index, slice):
+            index = _slice_slice(self._index[_vdim], allM, index, M)
+
+        if _vdim == 0:
+            row, col = index, self._index[1]
+        else:
+            row, col = self._index[0], index
+
+        if isinstance(row, INT_TYPES):
+            if isinstance(col, INT_TYPES):
+                return self._getitem_intXint(row, col)
+            elif isinstance(col, slice):
+                return self._getitem_intXslice(row, col)
+        elif isinstance(row, slice):
+            assert(isinstance(col, INT_TYPES))
+            return self._getitem_sliceXint(row, col)
+
+        raise IndexError("unsupported column index types")
+
+    def _getitem_intXint(self, row, col):
+        return self.m[row, col]
+
+    def _getitem_intXslice(self, row, col):
+        shape = (_slice_length(col, self.m.shape[1]), )
+        return self.__class__(self.m, shape=shape, index=(row, col))
+
+    def _getitem_sliceXint(self, row, col):
+        shape = (_slice_length(row, self.m.shape[0]), )
+        return self.__class__(self.m, shape=shape, index=(row, col))
+
+    def toarray(self):
+        return self.m[self._index]
+
+
+def _unpack_index(index, shape):
+    if not isinstance(index, tuple):
+        index = (index, )
+    if len(shape) < len(index):
+        raise IndexError("invalid index dimensionality - must be 2")
+
+    unpacked = ()
+    for shp, idx in zip_longest(shape, index):
+        idx = slice(None) if idx is None else idx
+        idx = _slice_defaults(idx, shp) if isinstance(idx, slice) else idx
+        unpacked += (idx, )
+
+    return unpacked
+
+
+def _slice_slice(outer, outer_len, inner, inner_len):
+    """
+    slice a slice - we take advantage of Python 3 range's support
+    for indexing.
+    """
+    assert(outer_len >= inner_len)
+    outer_rng = range(*outer.indices(outer_len))
+    rng = outer_rng[inner]
+    start, stop, step = rng.start, rng.stop, rng.step
+    if step < 0 and stop < 0:
+        stop = None
+    return slice(start, stop, step)
+
+
+def _range_length(start, stop, step):
+    """ return length of range """
+    assert(step != 0)
+    assert(start is not None and stop is not None and step is not None)
+    if step > 0 and start < stop:
+        return 1 + (stop - 1 - start) // step
+    elif step < 0 and start > stop:
+        return 1 + (start - 1 - stop) // -step
+    else:
+        return 0
+
+
+def _slice_length(s, length):
+    """ return slice length """
+    return _range_length(*s.indices(length))
+
+
+def _slice_defaults(s, length):
+    """ apply slice defaulting conventions """
+    assert(length >= 0)
+
+    step = 1 if s.step is None else s.step
+
+    if s.start is not None:
+        start = s.start
+        if start < 0:
+            start += length
+    else:
+        start = 0 if step > 0 else (length - 1)
+
+    if s.stop is not None:
+        stop = s.stop
+        if stop < 0:
+            stop += length
+    else:
+        stop = length if step > 0 else -length - 1
+
+    return slice(start, stop, step)

server/cli/launch.py

@@ -60,6 +60,20 @@ def common_args(func):
         metavar="<user labels CSV file>",
         help="CSV file containing user annotations; will be overwritten.  Created if does not exist.",
     )
+    @click.option(
+        "--backed",
+        is_flag=True,
+        default=False,
+        show_default=False,
+        help="Load data in file-backed mode, which may save memory, but result in slower overall performance."
+    )
+    @click.option(
+        "--disable-diffexp",
+        is_flag=True,
+        default=False,
+        show_default=False,
+        help="Disable on-demand differential expression."
+    )
     @functools.wraps(func)
     def wrapper(*args, **kwargs):
         return func(*args, **kwargs)
@@ -67,7 +81,9 @@ def common_args(func):
     return wrapper
 
 
-def parse_engine_args(embedding, obs_names, var_names, max_category_items, diffexp_lfc_cutoff, experimental_label_file):
+def parse_engine_args(embedding, obs_names, var_names, max_category_items,
+                      diffexp_lfc_cutoff, experimental_label_file, backed,
+                      disable_diffexp):
     return {
         "layout": embedding,
         "max_category_items": max_category_items,
@@ -75,6 +91,8 @@ def parse_engine_args(embedding, obs_names, var_names, max_category_items, diffe
         "obs_names": obs_names,
         "var_names": var_names,
         "label_file": experimental_label_file,
+        "backed": backed,
+        "disable_diffexp": disable_diffexp
     }
 
 
@@ -124,7 +142,9 @@ def launch(
         title,
         scripts,
         about,
-        experimental_label_file
+        experimental_label_file,
+        backed,
+        disable_diffexp
 ):
     """Launch the cellxgene data viewer.
     This web app lets you explore single-cell expression data.
@@ -140,7 +160,8 @@ def launch(
     > cellxgene launch <url>"""
 
     e_args = parse_engine_args(embedding, obs_names, var_names, max_category_items,
-                               diffexp_lfc_cutoff, experimental_label_file)
+                               diffexp_lfc_cutoff, experimental_label_file, backed,
+                               disable_diffexp)
     try:
         data_locator = DataLocator(data)
     except RuntimeError as re:
@@ -246,6 +267,10 @@ def launch(
     except ScanpyFileError as e:
         raise click.ClickException(f"{e}")
 
+    if not disable_diffexp and server.app.data.config['diffexp_may_be_slow']:
+        click.echo(f"[cellxgene] CAUTION: due to the size of your dataset, "
+                   f"running differential expression may take longer or fail.")
+
     if open_browser:
         click.echo(f"[cellxgene] Launching! Opening your browser to {cellxgene_url} now.")
         webbrowser.open(cellxgene_url)
@@ -259,7 +284,7 @@ def launch(
         sys.stdout = f
 
     try:
-        server.app.run(host=host, debug=debug, port=port, threaded=True, use_debugger=False)
+        server.app.run(host=host, debug=debug, port=port, threaded=False if debug else True, use_debugger=False)
     except OSError as e:
         if e.errno == errno.EADDRINUSE:
             raise click.ClickException("Port is in use, please specify an open port using the --port flag.") from e

server/requirements.txt

@@ -9,7 +9,7 @@ flatbuffers>=1.10.0
 fsspec>=0.4.4
 numpy>=1.15.2
 pandas>=0.23.1
-scipy>=1.1.0
+scipy>=1.3.0
 tables==3.5.1
 # TEMP workaround for https://github.com/theislab/scanpy/issues/832 aka h5py regression
 h5py==2.9.0

server/test/test_matrix_proxy.py

@@ -0,0 +1,233 @@
+import unittest
+import numpy as np
+from server.app.util.matrix_proxy import MatrixProxyView, MatrixProxy
+
+
+class NdArrayProxyView(MatrixProxyView):
+    """
+    Fake test class for matrix proxy - wraps ndarray
+    """
+    @classmethod
+    def __supports__(cls):
+        return ('numpy.ndarray', )
+
+
+class MatrixProxyViewTest(unittest.TestCase):
+    def test_ismatrixproxy(self):
+        n = np.zeros((2, 4))
+        mp = MatrixProxy.create(n)
+
+        self.assertIsNotNone(n)
+        self.assertIsNotNone(mp)
+        self.assertTrue(isinstance(mp, NdArrayProxyView))
+        self.assertFalse(MatrixProxy.ismatrixproxy(n))
+        self.assertTrue(MatrixProxy.ismatrixproxy(mp))
+
+    def test_params(self):
+        n = np.arange(15, dtype=np.float32).reshape((3, 5))
+        mp = MatrixProxy.create(n)
+
+        self.assertTrue(MatrixProxy.ismatrixproxy(mp))
+        self.assertEqual(mp.ndim, 2)
+        self.assertEqual(mp.shape, (3, 5))
+        self.assertEqual(mp.dtype, np.float32)
+
+        mpt = mp.T
+        self.assertTrue(MatrixProxy.ismatrixproxy(mpt))
+        self.assertEqual(mpt.ndim, 2)
+        self.assertEqual(mpt.shape, (5, 3))
+        self.assertEqual(mpt.dtype, np.float32)
+
+    def test_toarray(self):
+        n = np.arange(15, dtype=np.float32).reshape((3, 5))
+        mp = MatrixProxy.create(n)
+        self.assertTrue(np.all(mp.toarray() == [
+            [0., 1., 2., 3., 4.],
+            [5., 6., 7., 8., 9.],
+            [10., 11., 12., 13., 14.]
+        ]))
+        self.assertTrue(np.all(mp.T.toarray() == [
+            [0., 5., 10.],
+            [1., 6., 11.],
+            [2., 7., 12.],
+            [3., 8., 13.],
+            [4., 9., 14.]
+        ]))
+
+    def test_indexing(self):
+        """
+        [[ 0.,  1.,  2.,  3.,  4.],
+         [ 5.,  6.,  7.,  8.,  9.],
+         [10., 11., 12., 13., 14.]]
+        """
+        n = np.arange(15, dtype=np.float32).reshape((3, 5))
+        mp = MatrixProxy.create(n)
+
+        # should support int or slice for one or both dimensions
+
+        # int, int
+
+        self.assertEqual(mp[0, 0], 0)
+        self.assertEqual(mp.T[0, 0], 0)
+        self.assertEqual(mp[2, 4], 14)
+        self.assertEqual(mp.T[4, 2], 14)
+        self.assertEqual(mp[1, 3], mp.T[3, 1])
+
+        # int, slice
+
+        self.assertTrue(np.all(mp[0].toarray() == [0, 1, 2, 3, 4]))
+        self.assertTrue(np.all(mp[1, 1:].toarray() == [6, 7, 8, 9]))
+        self.assertTrue(np.all(mp[2, 1::-1].toarray() == [11, 10]))
+
+        self.assertTrue(np.all(mp.T[0].toarray() == [0, 5, 10]))
+        self.assertTrue(np.all(mp.T[1, 1:].toarray() == [6, 11]))
+        self.assertTrue(np.all(mp.T[2, 1::-1].toarray() == [7, 2]))
+
+        # slice, int
+
+        self.assertTrue(np.all(mp[:, 0].toarray() == [0, 5, 10]))
+        self.assertTrue(np.all(mp[1:, 1].toarray() == [6, 11]))
+        self.assertTrue(np.all(mp[1::-1, 2].toarray() == [7, 2]))
+
+        self.assertTrue(np.all(mp.T[:, 0].toarray() == [0, 1, 2, 3, 4]))
+        self.assertTrue(np.all(mp.T[1:, 1].toarray() == [6, 7, 8, 9]))
+        self.assertTrue(np.all(mp.T[1::-1, 2].toarray() == [11, 10]))
+
+        # slice, slice
+
+        self.assertTrue(np.all(mp[1:3, 2:4].toarray() == [
+            [7, 8],
+            [12, 13]
+        ]))
+        self.assertTrue(np.all(mp[:3, :4].toarray() == [
+            [0., 1., 2., 3.],
+            [5., 6., 7., 8.],
+            [10., 11., 12., 13.]
+        ]))
+        self.assertTrue(np.all(mp[::-1, ::-1].toarray() == [
+            [14, 13, 12, 11, 10],
+            [9, 8, 7, 6, 5],
+            [4, 3, 2, 1, 0]
+        ]))
+        self.assertTrue(np.all(mp[::-2, ::-2].toarray() == [
+            [14, 12, 10],
+            [4, 2, 0]
+        ]))
+
+        self.assertTrue(np.all(mp.T[2:4, 1:3].toarray() == [
+            [7, 12],
+            [8, 13]
+        ]))
+        self.assertTrue(np.all(mp.T[:4, :3].toarray() == [
+            [0, 5, 10],
+            [1, 6, 11],
+            [2, 7, 12],
+            [3, 8, 13]
+        ]))
+        self.assertTrue(np.all(mp.T[::-1, ::-1].toarray() == [
+            [14, 9, 4],
+            [13, 8, 3],
+            [12, 7, 2],
+            [11, 6, 1],
+            [10, 5, 0]
+        ]))
+        self.assertTrue(np.all(mp.T[::-2, ::-2].toarray() == [
+            [14, 4],
+            [12, 2],
+            [10, 0]
+        ]))
+
+    def test_repeated_indexing(self):
+        """
+        [[ 0.,  1.,  2.,  3.,  4.],
+         [ 5.,  6.,  7.,  8.,  9.],
+         [10., 11., 12., 13., 14.]]
+        """
+        n = np.arange(15, dtype=np.float32).reshape((3, 5))
+        mp = MatrixProxy.create(n)
+
+        self.assertEqual(mp[0][1], 1)
+        self.assertEqual(mp.T[0][1], 5)
+
+        self.assertTrue(np.all(mp[0::-1, ::-1][0, 2:4].toarray() == [
+            2, 1
+        ]))
+        self.assertTrue(np.all(mp[0::-1, 1:5:1][0, 1:3:1].toarray() == [
+            2, 3
+        ]))
+        self.assertTrue(np.all(mp[0::-1, 1:5:1][0, 2:0:-1].toarray() == [
+            3, 2
+        ]))
+        self.assertTrue(np.all(mp[0::-1, 5:1:-1][0, 1:3:1].toarray() == [
+            3, 2
+        ]))
+        self.assertTrue(np.all(mp[0::-1, 5:1:-1][0, 2:0:-1].toarray() == [
+            2, 3
+        ]))
+
+        self.assertTrue(np.all(mp.T[::-1, 0::-1][2:4, 0].toarray() == [
+            2, 1
+        ]))
+        self.assertTrue(np.all(mp.T[0::-1, 1:5:1][0, 1:3:1].toarray() == [
+            10
+        ]))
+        self.assertTrue(np.all(mp.T[0::-1, 1:5:1][0, 2:0:-1].toarray() == [
+            10
+        ]))
+        self.assertTrue(np.all(mp.T[0::-1, 5:1:-1][0, 1:3:1].toarray() == []))
+        self.assertTrue(np.all(mp.T[0::-1, 5:1:-1][0, 2:0:-1].toarray() == []))
+
+    def test_dimension_drop(self):
+        """
+        [[ 0.,  1.,  2.,  3.,  4.],
+         [ 5.,  6.,  7.,  8.,  9.],
+         [10., 11., 12., 13., 14.]]
+        """
+        n = np.arange(15, dtype=np.float32).reshape((3, 5))
+        mp = MatrixProxy.create(n)
+
+        # drop both dimensions, to a scalar
+        self.assertEqual(mp[0, 0], 0)
+
+        # drop 1 dimension, to an array
+        self.assertTrue(np.all(mp[0, :].toarray() == [
+            0, 1, 2, 3, 4
+        ]))
+        self.assertTrue(np.all(mp[:, 0].toarray() == [
+            0, 5, 10
+        ]))
+
+        # with .T
+        self.assertTrue(np.all(mp[0:2].T[-1:].toarray() == [
+            [4, 9]
+        ]))
+        self.assertTrue(np.all(mp[0:2].T[-1].toarray() == [
+            4, 9
+        ]))
+
+    def test_iter(self):
+        """
+        check that __iter__ is doing what we expect
+        """
+        n = np.arange(15, dtype=np.float32).reshape((3, 5))
+        mp = MatrixProxy.create(n)
+
+        rows = [r for r in mp]
+        self.assertEqual(len(rows), 3)
+        self.assertTrue(np.all(rows[0].toarray() == [
+            0, 1, 2, 3, 4
+        ]))
+        for i, r in enumerate(rows):
+            self.assertTrue(np.all(mp[i].toarray() == r.toarray()))
+
+        cols = [c for c in mp.T]
+        self.assertEqual(len(cols), 5)
+        self.assertTrue(np.all(cols[0].toarray() == [
+            0, 5, 10
+        ]))
+        for i, c in enumerate(cols):
+            self.assertTrue(np.all(mp.T[i].toarray() == c.toarray()))
+
+        e = [e for e in mp[0]]
+        self.assertEqual(len(e), 5)
+        self.assertEqual(e, [0, 1, 2, 3, 4])

server/test/test_scanpy_engine.py

@@ -18,10 +18,14 @@ Test the scanpy engine using the pbmc3k data set.
 """
 
 
-@parameterized_class(("data_locator",), [
+@parameterized_class(("data_locator", "backed"), [
-    ("example-dataset/pbmc3k.h5ad",),
+    ("example-dataset/pbmc3k.h5ad", False),
-    ("server/test/test_datasets/pbmc3k-CSC-gz.h5ad",),
+    ("server/test/test_datasets/pbmc3k-CSC-gz.h5ad", False),
-    ("server/test/test_datasets/pbmc3k-CSR-gz.h5ad",)
+    ("server/test/test_datasets/pbmc3k-CSR-gz.h5ad", False),
+
+    ("example-dataset/pbmc3k.h5ad", True),
+    ("server/test/test_datasets/pbmc3k-CSC-gz.h5ad", True),
+    ("server/test/test_datasets/pbmc3k-CSR-gz.h5ad", True),
 ])
 class EngineTest(unittest.TestCase):
     def setUp(self):
@@ -31,7 +35,8 @@ class EngineTest(unittest.TestCase):
             "obs_names": None,
             "var_names": None,
             "diffexp_lfc_cutoff": 0.01,
-            "layout_file": None
+            "layout_file": None,
+            "backed": self.backed
         }
         self.data = ScanpyEngine(DataLocator(self.data_locator), args)
 
@@ -51,9 +56,12 @@ class EngineTest(unittest.TestCase):
 
     @pytest.mark.filterwarnings("ignore:Scanpy data matrix")
     def test_data_type(self):
-        self.data.data.X = self.data.data.X.astype("float64")
+        # don't run the test on the more exotic data types, as they don't
-        with self.assertWarns(UserWarning):
+        # support the astype() interface (used by this test, but not underlying app)
-            self.data._validate_data_types()
+        if isinstance(self.data.data.X, np.ndarray):
+            self.data.data.X = self.data.data.X.astype("float64")
+            with self.assertWarns(UserWarning):
+                self.data._validate_data_types()
 
     def test_filter_idx(self):
         filter_ = {
@@ -159,10 +167,11 @@ class EngineTest(unittest.TestCase):
         self.assertEqual(len(result), 20)
 
     def test_data_frame(self):
-        fbs = self.data.data_frame_to_fbs_matrix(None, "var")
+        f1 = {"var": {"index": [[0, 10]]}}
+        fbs = self.data.data_frame_to_fbs_matrix(f1, "var")
         data = decode_fbs.decode_matrix_FBS(fbs)
         self.assertEqual(data["n_rows"], 2638)
-        self.assertEqual(data["n_cols"], 1838)
+        self.assertEqual(data["n_cols"], 10)
 
         with self.assertRaises(ValueError):
             self.data.data_frame_to_fbs_matrix(None, "obs")

server/test/test_scanpy_engine_data_load.py

@@ -25,6 +25,9 @@ class DataLoadEngineTest(unittest.TestCase):
             "var_names": "bar",
             "diffexp_lfc_cutoff": 0.1,
             "label_file": None,
+            "backed": False,
+            "diffexp_may_be_slow": False,
+            "disable_diffexp": False
         }
         self.data.update(args=args)
         self.assertEqual(args, self.data.config)