Merge branch 'master' into colinmegill/geneset-prototype

server/app/app.py

@@ -58,7 +58,7 @@ def cache_control_always(**cache_kwargs):
 
 @webbp.route("/", methods=["GET"])
 @cache_control(public=True, max_age=ONE_WEEK)
-def dataset_index(dataset=None):
+def dataset_index(url_dataroot=None, dataset=None):
     config = current_app.app_config
     if dataset is None:
         if config.single_dataset__datapath:
@@ -66,7 +66,10 @@ def dataset_index(dataset=None):
         else:
             return dataroot_index()
     else:
-        location = path_join(config.multi_dataset__dataroot, dataset)
+        dataroot = config.multi_dataset__dataroot.get(url_dataroot)
+        if dataroot is None:
+            abort(HTTPStatus.NOT_FOUND)
+        location = path_join(dataroot, dataset)
 
     scripts = config.server__scripts
     inline_scripts = config.server__inline_scripts
@@ -91,18 +94,21 @@ def health():
     return health_check(config)
 
 
-def get_data_adaptor(dataset=None):
+def get_data_adaptor(url_dataroot=None, dataset=None):
     config = current_app.app_config
 
     if dataset is None:
         datapath = config.single_dataset__datapath
     else:
-        datapath = path_join(config.multi_dataset__dataroot, dataset)
+        dataroot = config.multi_dataset__dataroot.get(url_dataroot)
+        if dataroot is None:
+            raise DatasetAccessError(f"Invalid dataset {url_dataroot}/{dataset}")
+        datapath = path_join(dataroot, dataset)
         # path_join returns a normalized path.  Therefore it is
         # sufficient to check that the datapath starts with the
         # dataroot to determine that the datapath is under the dataroot.
-        if not datapath.startswith(config.multi_dataset__dataroot):
-            raise DatasetAccessError("Invalid dataset {dataset}")
+        if not datapath.startswith(dataroot):
+            raise DatasetAccessError("Invalid dataset {url_dataroot}/{dataset}")
 
     if datapath is None:
         return common_rest.abort_and_log(HTTPStatus.BAD_REQUEST, "Invalid dataset NONE", loglevel=logging.INFO)
@@ -115,7 +121,7 @@ def rest_get_data_adaptor(func):
     @wraps(func)
     def wrapped_function(self, dataset=None):
         try:
-            with get_data_adaptor(dataset) as data_adaptor:
+            with get_data_adaptor(self.url_dataroot, dataset) as data_adaptor:
                 return func(self, data_adaptor)
         except DatasetAccessError:
             return common_rest.abort_and_log(
@@ -132,22 +138,23 @@ def dataroot_test_index():
     data += "<body><H1>Welcome to cellxgene</H1>"
 
     config = current_app.app_config
-    locator = DataLocator(config.multi_dataset__dataroot, region_name=config.data_locator__s3__region_name)
     datasets = []
-    for fname in locator.ls():
-        location = path_join(config.multi_dataset__dataroot, fname)
-        try:
-            MatrixDataLoader(location, app_config=config)
-            datasets.append(fname)
-        except DatasetAccessError:
-            # skip over invalid datasets
-            pass
+    for url_dataroot, dataroot in config.multi_dataset__dataroot.items():
+        locator = DataLocator(dataroot, region_name=config.data_locator__s3__region_name)
+        for fname in locator.ls():
+            location = path_join(dataroot, fname)
+            try:
+                MatrixDataLoader(location, app_config=config)
+                datasets.append((url_dataroot, fname))
+            except DatasetAccessError:
+                # skip over invalid datasets
+                pass
 
     data += "<br/>Select one of these datasets...<br/>"
     data += "<ul>"
     datasets.sort()
-    for dataset in datasets:
-        data += f"<li><a href=d/{dataset}>{dataset}</a></li>"
+    for url_dataroot, dataset in datasets:
+        data += f"<li><a href={url_dataroot}/{dataset}>{dataset}</a></li>"
     data += "</ul>"
     data += "</body></html>"
 
@@ -165,21 +172,29 @@ def dataroot_index():
         return redirect(config.multi_dataset__index)
 
 
-class SchemaAPI(Resource):
+class DatasetResource(Resource):
+    """Base class for all Resources that act on datasets."""
+
+    def __init__(self, url_dataroot):
+        super().__init__()
+        self.url_dataroot = url_dataroot
+
+
+class SchemaAPI(DatasetResource):
     @cache_control(public=True, max_age=ONE_WEEK)
     @rest_get_data_adaptor
     def get(self, data_adaptor):
         return common_rest.schema_get(data_adaptor, current_app.annotations)
 
 
-class ConfigAPI(Resource):
+class ConfigAPI(DatasetResource):
     @cache_control(public=True, max_age=ONE_WEEK)
     @rest_get_data_adaptor
     def get(self, data_adaptor):
         return common_rest.config_get(current_app.app_config, data_adaptor, current_app.annotations)
 
 
-class AnnotationsObsAPI(Resource):
+class AnnotationsObsAPI(DatasetResource):
     @cache_control(public=True, max_age=ONE_WEEK)
     @rest_get_data_adaptor
     def get(self, data_adaptor):
@@ -191,14 +206,14 @@ class AnnotationsObsAPI(Resource):
         return common_rest.annotations_obs_put(request, data_adaptor, current_app.annotations)
 
 
-class AnnotationsVarAPI(Resource):
+class AnnotationsVarAPI(DatasetResource):
     @cache_control(public=True, max_age=ONE_WEEK)
     @rest_get_data_adaptor
     def get(self, data_adaptor):
         return common_rest.annotations_var_get(request, data_adaptor, current_app.annotations)
 
 
-class DataVarAPI(Resource):
+class DataVarAPI(DatasetResource):
     @cache_control(no_store=True)
     @rest_get_data_adaptor
     def put(self, data_adaptor):
@@ -210,21 +225,21 @@ class DataVarAPI(Resource):
         return common_rest.data_var_get(request, data_adaptor)
 
 
-class ColorsAPI(Resource):
+class ColorsAPI(DatasetResource):
     @cache_control(public=True, max_age=ONE_WEEK)
     @rest_get_data_adaptor
     def get(self, data_adaptor):
         return common_rest.colors_get(data_adaptor)
 
 
-class DiffExpObsAPI(Resource):
+class DiffExpObsAPI(DatasetResource):
     @cache_control(no_store=True)
     @rest_get_data_adaptor
     def post(self, data_adaptor):
         return common_rest.diffexp_obs_post(request, data_adaptor)
 
 
-class LayoutObsAPI(Resource):
+class LayoutObsAPI(DatasetResource):
     @cache_control(public=True, max_age=ONE_WEEK)
     @rest_get_data_adaptor
     def get(self, data_adaptor):
@@ -236,20 +251,25 @@ class LayoutObsAPI(Resource):
         return common_rest.layout_obs_put(request, data_adaptor)
 
 
-def get_api_resources(bp_api):
+def get_api_resources(bp_api, url_dataroot=None):
     api = Api(bp_api)
+
+    def add_resource(resource, url):
+        """convenience function to make the outer function less verbose"""
+        api.add_resource(resource, url, resource_class_args=(url_dataroot,))
+
     # Initialization routes
-    api.add_resource(SchemaAPI, "/schema")
-    api.add_resource(ConfigAPI, "/config")
+    add_resource(SchemaAPI, "/schema")
+    add_resource(ConfigAPI, "/config")
     # Data routes
-    api.add_resource(AnnotationsObsAPI, "/annotations/obs")
-    api.add_resource(AnnotationsVarAPI, "/annotations/var")
-    api.add_resource(DataVarAPI, "/data/var")
+    add_resource(AnnotationsObsAPI, "/annotations/obs")
+    add_resource(AnnotationsVarAPI, "/annotations/var")
+    add_resource(DataVarAPI, "/data/var")
     # Display routes
-    api.add_resource(ColorsAPI, "/colors")
+    add_resource(ColorsAPI, "/colors")
     # Computation routes
-    api.add_resource(DiffExpObsAPI, "/diffexp/obs")
-    api.add_resource(LayoutObsAPI, "/layout/obs")
+    add_resource(DiffExpObsAPI, "/diffexp/obs")
+    add_resource(LayoutObsAPI, "/layout/obs")
     return api
 
 
@@ -285,10 +305,18 @@ class Server:
             # NOTE:  These routes only allow the dataset to be in the directory
             # of the dataroot, and not a subdirectory.  We may want to change
             # the route format at some point
-            bp_api = Blueprint("api_dataset", __name__, url_prefix="/d/<dataset>" + api_version)
-            resources = get_api_resources(bp_api)
-            self.app.register_blueprint(resources.blueprint)
-            self.app.add_url_rule("/d/<dataset>/", "dataset_index", dataset_index, methods=["GET"])
+            for url_dataroot in app_config.multi_dataset__dataroot.keys():
+                bp_api = Blueprint(
+                    f"api_dataset_{url_dataroot}", __name__, url_prefix=f"/{url_dataroot}/<dataset>" + api_version
+                )
+                resources = get_api_resources(bp_api, url_dataroot)
+                self.app.register_blueprint(resources.blueprint)
+                self.app.add_url_rule(
+                    f"/{url_dataroot}/<dataset>/",
+                    f"dataset_index_{url_dataroot}",
+                    lambda dataset: dataset_index(url_dataroot, dataset),
+                    methods=["GET"],
+                )
         self.app.matrix_data_cache_manager = app_config.matrix_data_cache_manager
         self.app.annotations = app_config.user_annotations
         self.app.app_config = app_config

server/common/app_config.py

@@ -1,9 +1,9 @@
 from server import __version__ as cellxgene_version
-from flatten_dict import flatten
+from flatten_dict import flatten, unflatten
 import os
 from os.path import splitext, basename, isdir
 import sys
-from urllib.parse import urlparse
+from urllib.parse import urlparse, quote_plus
 import yaml
 import copy
 
@@ -125,17 +125,23 @@ class AppConfig(object):
         dc = copy.deepcopy(config)
         mapping = {}
 
-        # special case for tiledb_ctx whose value is a dict.
-        val = config.get("adaptor", {}).get("cxg_adaptor", {}).get("tiledb_ctx")
-        if val is not None:
-            mapping["adaptor__cxg_adaptor__tiledb_ctx"] = (("adaptor", "cxg_adaptor", "tiledb_ctx"), val)
-            del dc["adaptor"]["cxg_adaptor"]["tiledb_ctx"]
-
-        # special case for csp_directives whose value is a dict.
-        val = config.get("server", {}).get("csp_directives")
-        if val is not None:
-            mapping["server__csp_directives"] = (("server", "csp_directives"), val)
-            del dc["server"]["csp_directives"]
+        # special cases where the value could be a dict.
+        # If its value is not None, the entry is added to the mapping, and not included
+        # in the flattening below.
+        dictval_cases = [
+            ("adaptor", "cxg_adaptor", "tiledb_ctx"),
+            ("server", "csp_directives"),
+            ("multi_dataset", "dataroot"),
+        ]
+        for dictval_case in dictval_cases:
+            cur = dc
+            for part in dictval_case[:-1]:
+                cur = cur.get(part, {})
+            val = cur.get(dictval_case[-1])
+            if val is not None:
+                key = "__".join(dictval_case)
+                mapping[key] = (dictval_case, val)
+                del cur[dictval_case[-1]]
 
         flat_config = flatten(dc)
         for key, value in flat_config.items():
@@ -162,6 +168,14 @@ class AppConfig(object):
 
         self.is_completed = False
 
+    def write_config(self, config_file):
+        """output the config to a yaml file"""
+        mapping = self.__mapping(self.default_config)
+        for attrname in mapping.keys():
+            mapping[attrname] = getattr(self, attrname)
+        config = unflatten(mapping, splitter=lambda key: key.split("__"))
+        yaml.dump(config, open(config_file, "w"))
+
     def update(self, **kw):
         for key, value in kw.items():
             if not hasattr(self, key):
@@ -302,6 +316,14 @@ class AppConfig(object):
         self.__check_attr("data_locator__s3__region_name", (type(None), bool, str))
         if self.data_locator__s3__region_name is True:
             path = self.single_dataset__datapath or self.multi_dataset__dataroot
+            if type(path) == dict:
+                # if multi_dataset__dataroot is a dict, then use the first key
+                # that is in s3.   NOTE:  it is not supported to have dataroots
+                # in different regions.
+                paths = path.values()
+                for path in paths:
+                    if path.startswith("s3://"):
+                        break
             if path.startswith("s3://"):
                 region_name = discover_s3_region_name(path)
                 if region_name is None:
@@ -366,7 +388,7 @@ class AppConfig(object):
                 )
 
     def handle_multi_dataset(self, context):
-        self.__check_attr("multi_dataset__dataroot", (type(None), str))
+        self.__check_attr("multi_dataset__dataroot", (type(None), dict, str))
         self.__check_attr("multi_dataset__index", (type(None), bool, str))
         self.__check_attr("multi_dataset__allowed_matrix_types", list)
         self.__check_attr("multi_dataset__matrix_cache__max_datasets", int)
@@ -375,6 +397,18 @@ class AppConfig(object):
         if self.multi_dataset__dataroot is None:
             return
 
+        if type(self.multi_dataset__dataroot) == str:
+            self.multi_dataset__dataroot = dict(d=self.multi_dataset__dataroot)
+
+        for key in self.multi_dataset__dataroot.keys():
+            # sanity check for well formed keys
+            if type(key) != str:
+                raise ConfigurationError(f"error in multi_dataset__dataroot {key}")
+            if quote_plus(key) != key:
+                raise ConfigurationError(f"error in multi_dataset__dataroot {key}")
+            if os.path.split(os.path.normpath(key))[-1] != key:
+                raise ConfigurationError(f"error in multi_dataset__dataroot {key}")
+
         # error checking
         for mtype in self.multi_dataset__allowed_matrix_types:
             try:

server/common/default_config.py

@@ -29,6 +29,23 @@ presentation:
   custom_colors: true
 
 multi_dataset:
+  # If dataroot is set, then cellxgene may serve multiple datasets.  This parameter is not
+  # compatable with single_dataset/datapath.
+  # dataroot may be a string, representing the path to a directory or S3 prefix.  In this
+  # case the datasets in that location are accessed from <server>/d/<datasetname>.
+  # example:
+  #     dataroot:  /path/to/datasets/
+  # or
+  #     dataroot:  s3://bucket/prefix/
+  #
+  # As an alternative, dataroot can be a dictionary, mapping url prefixes to dataroot paths.
+  # example:
+  #     dataroot:
+  #         set1 : /path/to/set1_datasets/
+  #         set2 : /path/to/set2_datasets/
+  # In this case, datasets can be accessed from <server>/set1/<datasetname> or
+  # <server>/set2/<datasetname>.
+
   dataroot: null
 
   # The index page when in multi-dataset mode:

server/common/health.py

@@ -23,12 +23,13 @@ def health_check(config):
     """
     health = {"status": None, "version": "1", "releaseID": cellxgene_version}
 
-    checks = [
-        (config.single_dataset__datapath is not None or config.multi_dataset__dataroot is not None),
-        _is_accessible(config.single_dataset__datapath, config),
-        _is_accessible(config.multi_dataset__dataroot, config),
-    ]
-    health["status"] = "pass" if all(checks) else "fail"
+    checks = False
+    if config.single_dataset__datapath is not None:
+        checks = _is_accessible(config.single_dataset__datapath, config)
+    elif config.multi_dataset__dataroot is not None:
+        checks = all([_is_accessible(datapath, config) for datapath in config.multi_dataset__dataroot.values()])
+
+    health["status"] = "pass" if checks else "fail"
     code = HTTPStatus.OK if health["status"] == "pass" else HTTPStatus.BAD_REQUEST
     response = make_response(jsonify(health), code)
     response.headers["Content-Type"] = "application/health+json"

server/compute/diffexp_cxg.py

@@ -3,6 +3,8 @@ import numpy as np
 from server.compute.diffexp_generic import diffexp_ttest_from_mean_var, mean_var_n
 from server.data_cxg.cxg_util import pack_selector_from_indices
 from server.common.errors import ComputeError
+from numba import jit
+
 
 """
 See the comments in diffexp_generic for a description of this algorithm
@@ -33,26 +35,33 @@ def get_thread_executor():
 
 
 def diffexp_ttest(adaptor, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
+
+    matrix = adaptor.open_array("X")
     row_selector_A = np.where(maskA)[0]
     row_selector_B = np.where(maskB)[0]
     nA = len(row_selector_A)
     nB = len(row_selector_B)
-    matrix = adaptor.open_array("X")
 
     dtype = matrix.dtype
     cols = matrix.shape[1]
     tile_extent = [dim.tile for dim in matrix.schema.domain]
 
-    # The rows from both row_selector_A and row_selector_B are gathered at the
-    # same time, then the mean and variance are computed by subsetting on that
-    # combined submatrix.  Combining the gather reduces number of requests/bandwidth
-    # to the data source.
-    row_selector_AB = np.union1d(row_selector_A, row_selector_B)
-    row_selector_A_in_AB = np.in1d(row_selector_AB, row_selector_A, assume_unique=True)
-    row_selector_B_in_AB = np.in1d(row_selector_AB, row_selector_B, assume_unique=True)
-    row_selector_AB = pack_selector_from_indices(row_selector_AB)
+    is_sparse = matrix.schema.sparse
 
-    # because all IO is done per-tile, and we are always dense and col-major,
+    if is_sparse:
+        row_selector_A = pack_selector_from_indices(row_selector_A)
+        row_selector_B = pack_selector_from_indices(row_selector_B)
+    else:
+        # The rows from both row_selector_A and row_selector_B are gathered at the
+        # same time, then the mean and variance are computed by subsetting on that
+        # combined submatrix.  Combining the gather reduces number of requests/bandwidth
+        # to the data source.
+        row_selector_AB = np.union1d(row_selector_A, row_selector_B)
+        row_selector_A_in_AB = np.in1d(row_selector_AB, row_selector_A, assume_unique=True)
+        row_selector_B_in_AB = np.in1d(row_selector_AB, row_selector_B, assume_unique=True)
+        row_selector_AB = pack_selector_from_indices(row_selector_AB)
+
+    # because all IO is done per-tile, and we are always col-major,
     # use the tile column size as the unit of partition.  Possibly access
     # more than one column tile at a time based on the target_workunit.
     # Revisit partitioning if we change the X layout, or start using a non-local execution environment
@@ -62,6 +71,7 @@ def diffexp_ttest(adaptor, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
     # the target_workunit.  A potential improvement would be to partition by both columns and rows.
     # However partitioning the rows is slightly more complex due to the arbitrary distribution
     # of row selections that are passed into this algorithm.
+
     cells_per_coltile = (nA + nB) * tile_extent[1]
     cols_per_partition = max(1, int(target_workunit / cells_per_coltile)) * tile_extent[1]
     col_partitions = [(c, min(c + cols_per_partition, cols)) for c in range(0, cols, cols_per_partition)]
@@ -73,10 +83,15 @@ def diffexp_ttest(adaptor, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
 
     executor = get_thread_executor()
     futures = []
-    for cols in col_partitions:
-        futures.append(
-            executor.submit(_mean_var_ab, matrix, row_selector_AB, row_selector_A_in_AB, row_selector_B_in_AB, cols)
-        )
+
+    if is_sparse:
+        for cols in col_partitions:
+            futures.append(executor.submit(_mean_var_sparse_ab, matrix, row_selector_A, nA, row_selector_B, nB, cols))
+    else:
+        for cols in col_partitions:
+            futures.append(
+                executor.submit(_mean_var_ab, matrix, row_selector_AB, row_selector_A_in_AB, row_selector_B_in_AB, cols)
+            )
 
     for future in futures:
         # returns tuple: (meanA, varA, meanB, varB, cols)
@@ -92,6 +107,12 @@ def diffexp_ttest(adaptor, maskA, maskB, top_n=8, diffexp_lfc_cutoff=0.01):
                 future.cancel()
             raise ComputeError(str(e))
 
+    if is_sparse:
+        if adaptor.has_array("X_col_shift"):
+            X_col_shift = adaptor.open_array("X_col_shift")[:]
+            meanA += X_col_shift
+            meanB += X_col_shift
+
     r = diffexp_ttest_from_mean_var(
         meanA.astype(dtype),
         varA.astype(dtype),
@@ -111,3 +132,67 @@ def _mean_var_ab(matrix, row_selector_AB, row_selector_A_in_AB, row_selector_B_i
     meanA, varA, n = mean_var_n(X[row_selector_A_in_AB])
     meanB, varB, n = mean_var_n(X[row_selector_B_in_AB])
     return (meanA, varA, meanB, varB, col_range)
+
+
+def _mean_var_sparse_ab(matrix, row_selector_A, nrows_A, row_selector_B, nrows_B, col_range):
+    meanA, varA = _mean_var_sparse(matrix, row_selector_A, nrows_A, col_range)
+    meanB, varB = _mean_var_sparse(matrix, row_selector_B, nrows_B, col_range)
+    return (meanA, varA, meanB, varB, col_range)
+
+
+@jit(nopython=True)
+def _mean_var_sparse_numba(x, var, nrows, ncols):
+    """Kernel to compute the mean and variance.  It was not clear if this function
+    could be written using numpy, thus avoiding the loops.  Therefore numba is
+    used here to speed things up.  With numba, this function takes a negligible amount
+    of time compared to reading in the sparse matrix"""
+    mean = np.zeros((ncols,), dtype=np.float64)
+    for col, val in zip(var, x):
+        mean[col] += val
+    mean /= nrows
+
+    # optimize the sumsq computation.
+    # since most entries in a sparse matrix are 0, then start by assuming
+    # all values are 0, so fill the sumsq array with nrows * (0 - mean)**2.
+    # as non-zero values are encountered, subtract off the (mean*mean) value
+    # and replace with (val-mean)**2.  Simplifying the expression
+    # gives the following code.
+    sumsq = nrows * np.multiply(mean, mean)
+    for col, val in zip(var, x):
+        sumsq[col] += val * (val - 2 * mean[col])
+    v = sumsq / (nrows - 1)
+    return mean, v
+
+
+def _mean_var_sparse(matrix, selector, nrows, col_range):
+    data = matrix.multi_index[selector, col_range[0] : col_range[1] - 1]
+    x = data[""]
+
+    # tiledb < 0.6.0 and >= 0.6.0 have slightly different interfaces.
+    # the following takes care of both cases:
+    #   older:  data["coords]["var"]
+    #   newer:  data["var"]
+    var = data.get("coords", data)["var"]
+
+    # shift the column indices to start at 0, this
+    # will become the index into the mean and var arrays.
+    var -= col_range[0]
+
+    fp_err_occurred = False
+
+    def fp_err_set(err, flag):
+        nonlocal fp_err_occurred
+        fp_err_occurred = True
+
+    ncols = col_range[1] - col_range[0]
+    with np.errstate(divide="call", invalid="call", call=fp_err_set):
+        mean, v = _mean_var_sparse_numba(x, var, nrows, ncols)
+
+    if fp_err_occurred:
+        mean[np.isfinite(mean) == False] = 0  # noqa: E712
+        v[np.isfinite(v) == False] = 0  # noqa: E712
+    else:
+        mean[np.isnan(mean)] = 0
+        v[np.isnan(v)] = 0
+
+    return mean, v

server/converters/cxgtool.py

@@ -8,8 +8,11 @@ The organization of the TileDB structure is:
     ├─ obs                          TileDB array containing cell (row) attributes, one attribute per
     │                               dataframe column, shape (n_obs,)
     ├─ var                          TileDB array containing gene (column) attributes, with one attribute per
-    │                               dataframe column, shape (n_obs,)
+    │                               dataframe column, shape (n_var,)
     ├─ X                            Main count matrix as a 2D TileDB array, single unnamed numeric attribute
+    ├─ X_col_shift                  TilebDB Array used in column shift encoding, shape (n_var,), dtype = X.dtype.
+    │                               Single unnamed numeric attribute.  If this array is sparse, and X_col_shift exists,
+    │                               then all values in the i'th column were subtracted by X_col_shift[i].
     ├─ emb                          TileDB group, storing optional embeddings (group may be empty)
     │  └─ <name1>                   TileDB Array, single anon attribute, ND numeric array, shape (n_obs, N)
     └─ cxg_group_metadata           Empty array used only to stash metadata about the overall object.
@@ -70,6 +73,7 @@ import argparse
 import numpy as np
 from os.path import splitext, basename
 import json
+from scipy.stats import mode
 
 from server.common.colors import convert_anndata_category_colors_to_cxg_category_colors
 from server.common.errors import ColorFormatException
@@ -114,6 +118,13 @@ def main():
         help="URL providing more information about the dataset (hint: must be a fully specified absolute URL).",
     )
     parser.add_argument("--out", "--output", "-o", help="output CXG file name")
+    parser.add_argument(
+        "--sparse-threshold",
+        "-s",
+        type=float,
+        default=0.0,  # force dense by default
+        help="The X array will be sparse if the percent of non-zeros falls below this value",
+    )
     args = parser.parse_args()
 
     global log_level
@@ -135,12 +146,15 @@ def main():
         obs_names=args.obs_names,
         about=args.about,
         extract_colors=not args.disable_custom_colors,
+        sparse_threshold=args.sparse_threshold,
     )
 
     log(1, "done")
 
 
-def write_cxg(adata, container, title, var_names=None, obs_names=None, about=None, extract_colors=False):
+def write_cxg(
+    adata, container, title, var_names=None, obs_names=None, about=None, extract_colors=False, sparse_threshold=5.0
+):
     if not adata.var.index.is_unique:
         raise ValueError("Variable index is not unique - unable to convert.")
     if not adata.obs.index.is_unique:
@@ -195,7 +209,7 @@ def write_cxg(adata, container, title, var_names=None, obs_names=None, about=Non
     log(1, "\t...embeddings created")
 
     # X matrix
-    save_X(container, adata, ctx)
+    save_X(container, adata.X, ctx, sparse_threshold)
     log(1, "\t...X created")
 
 
@@ -366,7 +380,7 @@ def create_emb(e_name, emb):
     dims = []
     for d in range(emb.ndim):
         shape = emb.shape
-        dims.append(tiledb.Dim("", domain=(0, shape[d] - 1), tile=min(shape[d], 1000), dtype=np.uint32))
+        dims.append(tiledb.Dim(domain=(0, shape[d] - 1), tile=min(shape[d], 1000), dtype=np.uint32))
     domain = tiledb.Domain(*dims)
     schema = tiledb.ArraySchema(
         domain=domain, sparse=False, attrs=attrs, capacity=1_000_000, cell_order="row-major", tile_order="row-major"
@@ -398,46 +412,175 @@ def save_embeddings(container, adata, ctx):
             log(1, f"\t\t...{name} embedding created")
 
 
-def create_X(X_name, shape):
+def create_X(X_name, shape, is_sparse):
     """
-    Dense, always.  Future task: explore if sparse encoding is worth the trouble
-    below a sparsity threshold.
-
-    The X matrix is access in both row and column oriented patterns, depending on the
+    The X matrix is accessed in both row and column oriented patterns, depending on the
     particular operation.  Because of the data type, default compression works best.
-    The tile size (50, 100) and global layout (row/col) was choosen empirically, by benchmarking
+    The tile size, (50, 100) for dense, and (512,2048) for sparse,
+    and global layout (row/col) was chosen empirically, by benchmarking
     the current cellxgene backend.
     """
     filters = tiledb.FilterList([tiledb.ZstdFilter()])
     attrs = [tiledb.Attr(dtype=np.float32, filters=filters)]
-    domain = tiledb.Domain(
-        tiledb.Dim(name="obs", domain=(0, shape[0] - 1), tile=min(shape[0], 50), dtype=np.uint32),
-        tiledb.Dim(name="var", domain=(0, shape[1] - 1), tile=min(shape[1], 100), dtype=np.uint32),
-    )
+    if is_sparse:
+        domain = tiledb.Domain(
+            tiledb.Dim(name="obs", domain=(0, shape[0] - 1), tile=min(shape[0], 512), dtype=np.uint32),
+            tiledb.Dim(name="var", domain=(0, shape[1] - 1), tile=min(shape[1], 2048), dtype=np.uint32),
+        )
+    else:
+        domain = tiledb.Domain(
+            tiledb.Dim(name="obs", domain=(0, shape[0] - 1), tile=min(shape[0], 50), dtype=np.uint32),
+            tiledb.Dim(name="var", domain=(0, shape[1] - 1), tile=min(shape[1], 100), dtype=np.uint32),
+        )
     schema = tiledb.ArraySchema(
-        domain=domain, sparse=False, attrs=attrs, cell_order="row-major", tile_order="col-major"
+        domain=domain, sparse=is_sparse, attrs=attrs, cell_order="row-major", tile_order="col-major"
     )
-    tiledb.DenseArray.create(X_name, schema)
+    if is_sparse:
+        tiledb.SparseArray.create(X_name, schema)
+    else:
+        tiledb.DenseArray.create(X_name, schema)
 
 
-def save_X(container, adata, ctx):
+def evaluate_for_sparse_encoding(xdata, sparse_threshold):
+    """
+    This function determines if the X matrix has a sparsity below the sparse_threshold.
+    This function also returns the number of non-zeros encountered and number
+    of elements evaluated.  This function may return before evaluating the whole X matrix
+    if it can be determined that X is not sparse enough.
+    """
+    shape = xdata.shape
+    stride = min(int(np.power(10, np.around(np.log10(1e9 / shape[1])))), 10_000)
+    nnz = 0
+    maxnnz = int(shape[0] * shape[1] * sparse_threshold / 100)
+    for row in range(0, shape[0], stride):
+        lim = min(row + stride, shape[0])
+        a = xdata[row:lim, :]
+        if type(a) is not np.ndarray:
+            a = a.toarray()
+        nnz += np.count_nonzero(a)
+        if nnz > maxnnz:
+            return (False, nnz, lim * shape[1])
+        log(2, "\t...rows", lim, "of", shape[0], "nnz", nnz, "nnz percent %5.2f%%" % (100 * nnz / (lim * shape[1])))
+
+    is_sparse = (100.0 * nnz / (shape[0] * shape[1])) < sparse_threshold
+    return (is_sparse, nnz, shape[0] * shape[1])
+
+
+def evaluate_for_sparse_column_shift_encoding(xdata, sparse_threshold):
+    """Column shift encoding works by taking the most common value in each column, then
+    subtracting that value from each element of the column.  If each column mostly contains
+    its most common value, then the resulting matrix can be very sparse.
+
+    This function determines if column shift encoding can be used to transform
+    the X matrix into a sparse matrix with a sparsity below the sparse_threshold.
+    If so, return the col_shift array that stores this encoding.
+    This function also returns the number of non-zeros encountered and number
+    of elements evaluated.  This function may return before evaluating the whole X matrix
+    if it can be determined that X cannot benefit from column shift encoding.
+    """
+    shape = xdata.shape
+    stride = max(1, 128_000_000 // shape[0])
+    col_shift = np.zeros(shape[1])
+    nnz = 0
+    maxnnz = int(shape[0] * shape[1] * sparse_threshold / 100)
+    for col in range(0, shape[1], stride):
+        lim = min(col + stride, shape[1])
+        a = xdata[:, col:lim]
+        if type(a) is not np.ndarray:
+            a = a.toarray()
+        m = mode(a)
+        col_shift[col:lim] = m.mode
+        nnz += shape[0] * (lim - col) - np.sum(m.count)
+        if nnz > maxnnz:
+            return (None, nnz, shape[0] * lim)
+        log(2, "\t...cols", lim, "of", shape[1], "nnz",
+            nnz, "nnz percent %5.2f%%" % (100 * nnz / (lim * shape[0])))
+
+    is_sparse = (100.0 * nnz / (shape[0] * shape[1])) < sparse_threshold
+    return (col_shift if is_sparse else None, nnz, shape[0] * shape[1])
+
+
+def save_X(container, xdata, ctx, sparse_threshold, expect_sparse=False):
     # Save X count matrix
     X_name = f"{container}/X"
-    shape = adata.X.shape
-    create_X(X_name, shape)
 
+    shape = xdata.shape
+    log(1, "\t...shape:", str(shape))
+
+    col_shift = None
+    if sparse_threshold == 100:
+        is_sparse = True
+    elif sparse_threshold == 0:
+        is_sparse = False
+    else:
+        is_sparse, nnz, nelem = evaluate_for_sparse_encoding(xdata, sparse_threshold)
+        percent = 100.0 * nnz / nelem
+        if nelem != shape[0] * shape[1]:
+            log(1, "\t...sparse=", is_sparse, "non-zeros percent (estimate): %6.2f" % percent)
+        else:
+            log(1, "\t...sparse=", is_sparse, "non-zeros:", nnz, "percent: %6.2f" % percent)
+
+        is_sparse = percent < sparse_threshold
+        if not is_sparse:
+            col_shift, nnz, nelem = evaluate_for_sparse_column_shift_encoding(xdata, sparse_threshold)
+            is_sparse = col_shift is not None
+            percent = 100.0 * nnz / nelem
+            if nelem != shape[0] * shape[1]:
+                log(1, "\t...sparse=", is_sparse, "col shift non-zeros percent (estimate): %6.2f" % percent)
+            else:
+                log(1, "\t...sparse=", is_sparse, "col shift non-zeros:", nnz, "percent: %6.2f" % percent)
+
+    if expect_sparse is True and is_sparse is False:
+        return False
+
+    create_X(X_name, shape, is_sparse)
     stride = min(int(np.power(10, np.around(np.log10(1e9 / shape[1])))), 10_000)
-    with tiledb.DenseArray(X_name, mode="w", ctx=ctx) as X:
-        for row in range(0, shape[0], stride):
-            lim = min(row + stride, shape[0])
-            a = adata.X[row:lim, :]
-            if type(a) is not np.ndarray:
-                a = a.toarray()
-            X[row:lim, :] = a
-            log(2, "\t...rows", row, "to", lim)
-        tiledb.consolidate(X_name, ctx=ctx)
+    if is_sparse:
+        if col_shift is not None:
+            log(1, "\t...output X as sparse matrix with column shift encoding")
+            X_col_shift_name = f"{container}/X_col_shift"
+            filters = tiledb.FilterList([tiledb.ZstdFilter()])
+            attrs = [tiledb.Attr(dtype=np.float32, filters=filters)]
+            domain = tiledb.Domain(tiledb.Dim(domain=(0, shape[1] - 1), tile=min(shape[1], 5000), dtype=np.uint32))
+            schema = tiledb.ArraySchema(domain=domain, attrs=attrs)
+            tiledb.DenseArray.create(X_col_shift_name, schema)
+            with tiledb.DenseArray(X_col_shift_name, mode="w", ctx=ctx) as X_col_shift:
+                X_col_shift[:] = col_shift
+            tiledb.consolidate(X_col_shift_name, ctx=ctx)
+        else:
+            log(1, "\t...output X as sparse matrix")
+
+        with tiledb.SparseArray(X_name, mode="w", ctx=ctx) as X:
+            nnz = 0
+            for row in range(0, shape[0], stride):
+                lim = min(row + stride, shape[0])
+                a = xdata[row:lim, :]
+                if type(a) is not np.ndarray:
+                    a = a.toarray()
+                if col_shift is not None:
+                    a = a - col_shift
+                indices = np.nonzero(a)
+                trow = indices[0] + row
+                nnz += indices[0].shape[0]
+                X[trow, indices[1]] = a[indices[0], indices[1]]
+                log(2, "\t...rows", lim, "of", shape[0], "nnz", nnz, "sparse", nnz / (lim * shape[1]))
+
+    else:
+        log(1, "\t...output X as dense matrix")
+        with tiledb.DenseArray(X_name, mode="w", ctx=ctx) as X:
+            for row in range(0, shape[0], stride):
+                lim = min(row + stride, shape[0])
+                a = xdata[row:lim, :]
+                if type(a) is not np.ndarray:
+                    a = a.toarray()
+                X[row:lim, :] = a
+                log(2, "\t...rows", row, "to", lim)
 
     tiledb.consolidate(X_name, ctx=ctx)
+    if hasattr(tiledb, "vacuum"):
+        tiledb.vacuum(X_name)
+
+    return is_sparse
 
 
 def save_metadata(container, metadata_dict):

server/converters/to_sparse.py

@@ -0,0 +1,62 @@
+"""
+Script to create a sparse dataset in CXG format based on an input dataset in CXG format.
+The input dataset is not modified.
+"""
+import os
+import shutil
+import tiledb
+import argparse
+import sys
+import server.converters.cxgtool as cxgtool
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("input", help="input cxg directory")
+    parser.add_argument("output", help="output cxg directory")
+    parser.add_argument("--overwrite", action="store_true", help="replace output cxg directory")
+    parser.add_argument("--verbose", "-v", action="count", default=0, help="verbose output")
+    parser.add_argument(
+        "--sparse-threshold",
+        "-s",
+        type=float,
+        default=5.0,  # default is 5% non-zero values
+        help="The X array will be sparse if the percent of non-zeros falls below this value",
+    )
+    args = parser.parse_args()
+
+    if os.path.exists(args.output):
+        print("output dir exists:", args.output)
+        if args.overwrite:
+            print("output dir removed:", args.output)
+            shutil.rmtree(args.output)
+        else:
+            print("use the overwrite option to remove the output directory")
+            sys.exit(1)
+
+    if not os.path.isdir(args.input):
+        print("input is not a directory", args.input)
+        sys.exit(1)
+
+    shutil.copytree(args.input, args.output,
+                    ignore=shutil.ignore_patterns("X", "X_col_shift"))
+
+    ctx = tiledb.Ctx(
+        {
+            "sm.num_reader_threads": 32,
+            "sm.num_writer_threads": 32,
+            "sm.consolidation.buffer_size": 1 * 1024 * 1024 * 1024,
+        }
+    )
+
+    with tiledb.DenseArray(os.path.join(args.input, "X"), mode="r", ctx=ctx) as X_in:
+        is_sparse = cxgtool.save_X(args.output, X_in, ctx, args.sparse_threshold, expect_sparse=True)
+
+    if is_sparse is False:
+        print("The array is not sparse, cleaning up, abort.")
+        shutil.rmtree(args.output)
+        sys.exit(1)
+
+
+if __name__ == "__main__":
+    main()

server/data_common/data_adaptor.py

@@ -241,8 +241,7 @@ class DataAdaptor(metaclass=ABCMeta):
         duplicate_columns = list(set(labels_df.columns) & set(obs_columns))
         if len(duplicate_columns) > 0:
             raise KeyError(
-                "Labels file may not contain column names which overlap "
-                f"with h5ad obs columns {duplicate_columns}"
+                "Labels file may not contain column names which overlap " f"with h5ad obs columns {duplicate_columns}"
             )
 
         # labels must have same count as obs annotations

server/data_cxg/cxg_adaptor.py

@@ -133,6 +133,10 @@ class CxgAdaptor(DataAdaptor):
             return False
         return True
 
+    def has_array(self, name):
+        a_type = tiledb.object_type(path_join(self.url, name), ctx=self.tiledb_ctx)
+        return a_type == "array"
+
     def _validate_and_initialize(self):
         """
         remember, preload_validation() has already been called, so
@@ -147,13 +151,7 @@ class CxgAdaptor(DataAdaptor):
         * version 0.1 -- metadata attache to cxg_group_metadata array.
           Same as 0, except it adds group metadata.
         """
-        a_type = tiledb.object_type(path_join(self.url, "cxg_group_metadata"), ctx=self.tiledb_ctx)
-        if a_type is None:
-            # version 0
-            cxg_version = "0.0"
-            title = None
-            about = None
-        elif a_type == "array":
+        if self.has_array("cxg_group_metadata"):
             # version >0
             gmd = self.open_array("cxg_group_metadata")
             cxg_version = gmd.meta["cxg_version"]
@@ -161,6 +159,11 @@ class CxgAdaptor(DataAdaptor):
                 cxg_properties = json.loads(gmd.meta["cxg_properties"])
                 title = cxg_properties.get("title", None)
                 about = cxg_properties.get("about", None)
+        else:
+            # version 0
+            cxg_version = "0.0"
+            title = None
+            about = None
 
         if cxg_version not in ["0.0", "0.1"]:
             raise DatasetAccessError(f"cxg matrix is not valid: {self.url}")
@@ -171,7 +174,11 @@ class CxgAdaptor(DataAdaptor):
 
     @staticmethod
     def _open_array(uri, tiledb_ctx):
-        return tiledb.DenseArray(uri, mode="r", ctx=tiledb_ctx)
+        with tiledb.Array(uri, mode="r", ctx=tiledb_ctx) as array:
+            if array.schema.sparse:
+                return tiledb.SparseArray(uri, mode="r", ctx=tiledb_ctx)
+            else:
+                return tiledb.DenseArray(uri, mode="r", ctx=tiledb_ctx)
 
     def open_array(self, name):
         try:
@@ -200,15 +207,73 @@ class CxgAdaptor(DataAdaptor):
         meta = self.open_array("cxg_group_metadata").meta
         return json.loads(meta["cxg_category_colors"]) if "cxg_category_colors" in meta else dict()
 
+    def __remap_indices(self, coord_range, coord_mask, coord_data):
+        """
+        This function maps the indices in coord_data, which could be in the range [0,coord_range), to
+        a range that only includes the number of indices encoded in coord_mask.
+        coord_range is the maxinum size of the range (e.g. get_shape()[0] or get_shape()[1])
+        coord_mask is a mask passed into the get_X_array, of size coord_range
+        coord_data are indices representing locations of non-zero values, in the range [0,coord_range).
+
+        For example, say
+        coord_mask = [1,0,1,0,0,1]
+        coord_data = [2,0,2,2,5]
+
+        The function computes the following:
+        indices = [0,2,5]
+        ncoord = 3
+        maprange = [0,1,2]
+        mapindex = [0,0,1,0,0,2]
+        coordindices = [1,0,1,1,2]
+        """
+        if coord_mask is None:
+            return coord_range, coord_data
+
+        indices = np.where(coord_mask)[0]
+        ncoord = indices.shape[0]
+        maprange = np.arange(ncoord)
+        mapindex = np.zeros(indices[-1] + 1, dtype=int)
+        mapindex[indices] = maprange
+        coordindices = mapindex[coord_data]
+        return ncoord, coordindices
+
     def get_X_array(self, obs_mask=None, var_mask=None):
         obs_items = pack_selector_from_mask(obs_mask)
         var_items = pack_selector_from_mask(var_mask)
+        if obs_items is None or var_items is None:
+            # If either zero rows or zero columns were selected, return an empty 2d array.
+            shape = self.get_shape()
+            obs_size = 0 if obs_items is None else shape[0] if obs_mask is None else np.count_nonzero(obs_mask)
+            var_size = 0 if var_items is None else shape[1] if var_mask is None else np.count_nonzero(var_mask)
+            return np.ndarray((obs_size, var_size))
+
         X = self.open_array("X")
-        if obs_items == slice(None) and var_items == slice(None):
-            data = X[:, :]
+
+        if X.schema.sparse:
+            if obs_items == slice(None) and var_items == slice(None):
+                data = X[:, :]
+            else:
+                data = X.multi_index[obs_items, var_items]
+
+            nrows, obsindices = self.__remap_indices(X.shape[0], obs_mask, data.get("coords", data)["obs"])
+            ncols, varindices = self.__remap_indices(X.shape[1], var_mask, data.get("coords", data)["var"])
+            densedata = np.zeros((nrows, ncols), dtype=self.get_X_array_dtype())
+            densedata[obsindices, varindices] = data[""]
+            if self.has_array("X_col_shift"):
+                X_col_shift = self.open_array("X_col_shift")
+                if var_items == slice(None):
+                    densedata += X_col_shift[:]
+                else:
+                    densedata += X_col_shift.multi_index[var_items][""]
+
+            return densedata
+
         else:
-            data = X.multi_index[obs_items, var_items][""]
-        return data
+            if obs_items == slice(None) and var_items == slice(None):
+                data = X[:, :]
+            else:
+                data = X.multi_index[obs_items, var_items][""]
+            return data
 
     def get_shape(self):
         X = self.open_array("X")

server/data_cxg/cxg_util.py

@@ -20,7 +20,7 @@ def pack_selector_from_mask(boolarray):
 def pack_selector_from_indices(selector):
 
     if len(selector) == 0:
-        return slice(None)
+        return None
 
     result = []
     current = slice(selector[0], selector[0])

server/eb/.gitignore

@@ -0,0 +1,5 @@
+
+# Elastic Beanstalk Files
+.elasticbeanstalk/*
+!.elasticbeanstalk/*.cfg.yml
+!.elasticbeanstalk/*.global.yml

server/eb/app.py

@@ -67,7 +67,6 @@ class WSGIServer(Server):
             "object-src": ["'none'"],
             "base-uri": ["'none'"],
             "frame-ancestors": ["'none'"],
-            "require-trusted-types-for": ["'script'"],
         }
 
         if not app.debug:

server/requirements.txt

@@ -11,12 +11,13 @@ flask-talisman>=0.7.0
 flatbuffers>=1.10.0
 flatten-dict>=0.2.0
 fsspec>=0.4.4
+numba>=0.49.1
 numpy>=1.16.0
 packaging>=20.0
 pandas>=0.24.2
 PyYAML>=5.3
 scipy>=1.3.0
 requests>=2.22.0
-tiledb==0.5.9
+tiledb>=0.5.9,!=0.6.0
 s3fs>=0.4.2
 gunicorn>=20.0.4

server/test/__init__.py

@@ -2,13 +2,19 @@ import random
 import shutil
 import string
 import tempfile
+import requests
+import time
+import os
+from subprocess import Popen
 from os import path, popen
+from contextlib import contextmanager
 
 import pandas as pd
 
 from server.common.annotations import AnnotationsLocalFile
 from server.common.data_locator import DataLocator
-from server.common.app_config import AppConfig
+from server.common.app_config import AppConfig, DEFAULT_SERVER_PORT
+from server.common.utils import find_available_port
 from server.data_common.fbs.matrix import encode_matrix_fbs
 from server.data_common.matrix_loader import MatrixDataLoader, MatrixDataType
 
@@ -60,7 +66,7 @@ def skip_if(condition, reason: str):
     return decorator
 
 
-def app_config(data_locator, backed=False):
+def app_config(data_locator, backed=False, extra={}):
     args = {
         "embeddings__names": ["umap", "tsne", "pca"],
         "presentation__max_categories": 100,
@@ -74,9 +80,61 @@ def app_config(data_locator, backed=False):
     }
     config = AppConfig()
     config.update(**args)
+    config.update(**extra)
     config.complete_config()
     return config
 
 
 def random_string(n):
     return "".join(random.choice(string.ascii_letters) for _ in range(n))
+
+
+@contextmanager
+def test_server(command_line_args=[], app_config=None):
+    """A context to run the cellxgene server.
+    Command line arguments can be passed in, as well as an app_config.
+    This function is meant to be used like this, for example:
+
+    with test_server(...) as server:
+       r = requests.get(f"{server}/...")
+       // check r
+
+    where the server can be accessed within the context, and is terminated when
+    the context is exited.
+    The port is automatically set using find_available_port.
+    The verbose flag is automatically set to True.
+    If an app_config is provided, then this function writes a temporary
+    yaml config file, which this server will read and parse.
+    """
+
+    port = DEFAULT_SERVER_PORT
+    port = find_available_port("localhost", port)
+    command = ["cellxgene", "--no-upgrade-check", "launch", "--verbose", "--port=%d" % port] + command_line_args
+
+    tempdir = None
+    if app_config:
+        tempdir = tempfile.TemporaryDirectory()
+        config_file = os.path.join(tempdir.name, "config.yaml")
+        app_config.write_config(config_file)
+        command.extend(["-c", config_file])
+
+    server = f"http://localhost:{port}"
+    ps = Popen(command)
+
+    for _ in range(10):
+        try:
+            requests.get(f"{server}/health")
+            break
+        except requests.exceptions.ConnectionError:
+            time.sleep(1)
+
+    if tempdir:
+        tempdir.cleanup()
+
+    try:
+        yield server
+    finally:
+        try:
+            ps.terminate()
+        except ProcessLookupError:
+            pass

server/test/create_test_matrix.py

@@ -0,0 +1,44 @@
+import anndata
+import argparse
+import random
+import scipy
+import numpy as np
+
+
+def main():
+    parser = argparse.ArgumentParser("A command to generate test h5ad files")
+    parser.add_argument("output", help="Name of the output file")
+    parser.add_argument("nobs", type=int, help="Number of observations (rows)")
+    parser.add_argument("nvar", type=int, help="Number of variables (columns)")
+    parser.add_argument("-n", "--nnz-percent", type=float, default=100, help="percent of non-zeros")
+    parser.add_argument("-c", "--col-shift", action="store_true", help="add a random value to each column")
+    parser.add_argument("--seed", type=int, default=None, help="add a random value to each column")
+
+    args = parser.parse_args()
+    create_test_h5ad(args.output, args.nobs, args.nvar, args.nnz_percent, args.col_shift, args.seed)
+
+
+def create_test_h5ad(outfile, nobs, nvar, nnz_percent=100, apply_col_shift=False, seed=None):
+    random.seed(seed)
+    np.random.seed(seed)
+    x = create_X_array(nobs, nvar, nnz_percent, apply_col_shift)
+    obsm = {"X_random": np.random.rand(nobs, 2).astype(np.float32)}
+    adata = anndata.AnnData(x, obsm=obsm)
+    adata.write(outfile)
+
+
+def create_X_array(nobs, nvar, nnz_percent, apply_col_shift):
+    if nnz_percent < 100:
+        array = scipy.sparse.random(nobs, nvar, nnz_percent * 0.01, dtype=np.float32, format="csc")
+    else:
+        array = np.random.rand(nobs, nvar).astype(np.float32)
+
+    if apply_col_shift:
+        col_shift = np.random.rand((nvar))
+        array += col_shift
+
+    return array
+
+
+if __name__ == "__main__":
+    main()

server/test/run_diffexp.py

@@ -15,8 +15,10 @@ from server.data_cxg.cxg_adaptor import CxgAdaptor
 def main():
     parser = argparse.ArgumentParser("A command to test diffexp")
     parser.add_argument("dataset", help="name of a dataset to load")
-    parser.add_argument("-na", "--numA", type=int, required=True, help="number of rows in group A")
-    parser.add_argument("-nb", "--numB", type=int, required=True, help="number of rows in group B")
+    parser.add_argument("-na", "--numA", type=int, help="number of rows in group A")
+    parser.add_argument("-nb", "--numB", type=int, help="number of rows in group B")
+    parser.add_argument("-va", "--varA", help="obs variable:value to use for group A")
+    parser.add_argument("-vb", "--varB", help="obs variable:value to use for group B")
     parser.add_argument("-t", "--trials", default=1, type=int, help="number of trials")
     parser.add_argument(
         "-a", "--alg", choices=("default", "generic", "cxg"), default="default", help="algorithm to use"
@@ -41,22 +43,34 @@ def main():
         if isinstance(adaptor, CxgAdaptor):
             adaptor.open_array("X").schema.dump()
 
-    numA = args.numA
-    numB = args.numB
+    random.seed(args.seed)
+    np.random.seed(args.seed)
     rows = adaptor.get_shape()[0]
 
-    random.seed(args.seed)
+    if args.numA:
+        filterA = random.sample(range(rows), args.numA)
+    elif args.varA:
+        vname, vval = args.varA.split(":")
+        filterA = get_filter_from_obs(adaptor, vname, vval)
+    else:
+        print("must supply numA or varA")
+        sys.exit(1)
 
-    if not args.new_selection:
-        samples = random.sample(range(rows), numA + numB)
-        filterA = samples[:numA]
-        filterB = samples[numA:]
+    if args.numB:
+        filterB = random.sample(range(rows), args.numB)
+    elif args.varB:
+        vname, vval = args.varB.split(":")
+        filterB = get_filter_from_obs(adaptor, vname, vval)
+    else:
+        print("must supply numB or varB")
+        sys.exit(1)
 
     for i in range(args.trials):
         if args.new_selection:
-            samples = random.sample(range(rows), numA + numB)
-            filterA = samples[:numA]
-            filterB = samples[numA:]
+            if args.numA:
+                filterA = random.sample(range(rows), args.numA)
+            if args.numB:
+                filterB = random.sample(range(rows), args.numB)
 
         maskA = np.zeros(rows, dtype=bool)
         maskA[filterA] = True
@@ -82,5 +96,22 @@ def main():
             print(res)
 
 
+def get_filter_from_obs(adaptor, obsname, obsval):
+    attrs = adaptor.get_obs_columns()
+    if obsname not in attrs:
+        print(f"Unknown obs attr {obsname}: expected on of {attrs}")
+        sys.exit(1)
+    obsvals = adaptor.query_obs_array(obsname)[:]
+    obsval = type(obsvals[0])(obsval)
+
+    vfilter = np.where(obsvals == obsval)[0]
+    if len(vfilter) == 0:
+        u = np.unique(obsvals)
+        print(f"Unknown value in variable {obsname}:{obsval}: expected one of {list(u)}")
+        sys.exit(1)
+
+    return vfilter
+
+
 if __name__ == "__main__":
     main()

server/test/test_api.py

@@ -244,6 +244,19 @@ class EndPoints(object):
         self.assertEqual(df["n_rows"], 2638)
         self.assertEqual(df["n_cols"], 1)
 
+    def test_data_get_unknown_filter_fbs(self):
+        index_col_name = self.schema["schema"]["annotations"]["var"]["index"]
+        endpoint = "data/var"
+        query = f"var:{index_col_name}=UNKNOWN"
+        url = f"{self.URL_BASE}{endpoint}?{query}"
+        header = {"Accept": "application/octet-stream"}
+        result = self.session.get(url, headers=header)
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+        self.assertEqual(result.headers["Content-Type"], "application/octet-stream")
+        df = decode_fbs.decode_matrix_FBS(result.content)
+        self.assertEqual(df["n_rows"], 2638)
+        self.assertEqual(df["n_cols"], 0)
+
     def test_data_put_single_var(self):
         endpoint = "data/var"
         url = f"{self.URL_BASE}{endpoint}"

server/test/test_app_config.py

@@ -1,5 +1,8 @@
 import unittest
 from server.common.app_config import AppConfig
+from server.common.errors import ConfigurationError
+from server.test import PROJECT_ROOT, test_server
+import requests
 
 # NOTE, there are more tests that should be written for AppConfig.
 # this is just a start.
@@ -26,3 +29,43 @@ class AppConfigTest(unittest.TestCase):
         c.update(server__scripts=("a", "b"), server__inline_scripts=["c", "d"])
         v = c.changes_from_default()
         self.assertCountEqual(v, [("server__scripts", ["a", "b"], []), ("server__inline_scripts", ["c", "d"], [])])
+
+    def test_multi_dataset(self):
+
+        c = AppConfig()
+        # test for illegal url_dataroots
+        for illegal in ("a/b", "../b", "!$*", "\\n", "", "(bad)"):
+            c.update(multi_dataset__dataroot={illegal: f"{PROJECT_ROOT}/example-dataset"})
+            with self.assertRaises(ConfigurationError):
+                c.complete_config()
+
+        # test for legal url_dataroots
+        for legal in (
+            "d",
+            "this.is-okay_",
+        ):
+            c.update(multi_dataset__dataroot={legal: f"{PROJECT_ROOT}/example-dataset"})
+            c.complete_config()
+
+        # test that multi dataroots work end to end
+        c.update(
+            multi_dataset__dataroot=dict(
+                set1=f"{PROJECT_ROOT}/example-dataset",
+                set2=f"{PROJECT_ROOT}/server/test/test_datasets"
+            )
+        )
+        c.complete_config()
+
+        with test_server(app_config=c) as server:
+            session = requests.Session()
+
+            r = session.get(f"{server}/set1/pbmc3k.h5ad/api/v0.2/config")
+            data_config = r.json()
+            assert data_config["config"]["displayNames"]["dataset"] == "pbmc3k"
+
+            r = session.get(f"{server}/set2/pbmc3k.cxg/api/v0.2/config")
+            data_config = r.json()
+            assert data_config["config"]["displayNames"]["dataset"] == "pbmc3k"
+
+            r = session.get(f"{server}/health")
+            assert r.json()["status"] == "pass"

server/test/test_diffexp.py

@@ -1,24 +1,24 @@
 import unittest
 from server.data_common.matrix_loader import MatrixDataLoader
-from server.common.app_config import AppConfig
+from server.test import PROJECT_ROOT, app_config
 import server.compute.diffexp_cxg as diffexp_cxg
 import server.compute.diffexp_generic as diffexp_generic
+from server.converters.cxgtool import write_cxg
+from server.test.create_test_matrix import create_test_h5ad
+from server.data_common.fbs.matrix import encode_matrix_fbs, decode_matrix_fbs
 import numpy as np
-
-from server.test import PROJECT_ROOT
+import tempfile
+import os
 
 
 class DiffExpTest(unittest.TestCase):
     """Tests the diffexp returns the expected results for one test case, using different
     adaptor types and different algorithms."""
 
-    def load_dataset(self, path):
-        app_config = AppConfig()
-        app_config.single_dataset__datapath = path
-        app_config.server__verbose = True
-        app_config.complete_config()
+    def load_dataset(self, path, extra={}):
+        config = app_config(path, extra=extra)
         loader = MatrixDataLoader(path)
-        adaptor = loader.open(app_config)
+        adaptor = loader.open(config)
         return adaptor
 
     def get_mask(self, adaptor, start, stride):
@@ -29,6 +29,14 @@ class DiffExpTest(unittest.TestCase):
         mask[sel] = True
         return mask
 
+    def compare_diffexp_results(self, results, expects):
+        self.assertEqual(len(results), len(expects))
+        for result, expect in zip(results, expects):
+            self.assertEqual(result[0], expect[0])
+            self.assertTrue(np.isclose(result[1], expect[1], 1e-6, 1e-4))
+            self.assertTrue(np.isclose(result[2], expect[2], 1e-6, 1e-4))
+            self.assertTrue(np.isclose(result[3], expect[3], 1e-6, 1e-4))
+
     def check_1_10_2_10(self, results):
         """Checks the results for a specific set of rows selections"""
         expects = [
@@ -43,12 +51,12 @@ class DiffExpTest(unittest.TestCase):
             [1575, 1.0317602, 0.007830310753043345, 1.0],
             [576, 0.97873515, 0.008272092578813124, 1.0],
         ]
-        self.assertEqual(len(results), len(expects))
-        for result, expect in zip(results, expects):
-            self.assertEqual(result[0], expect[0])
-            self.assertAlmostEqual(result[1], expect[1])
-            self.assertAlmostEqual(result[2], expect[2])
-            self.assertAlmostEqual(result[3], expect[3])
+        self.compare_diffexp_results(results, expects)
+
+    def get_X_col(self, adaptor, cols):
+        varmask = np.zeros(adaptor.get_shape()[1], dtype=bool)
+        varmask[cols] = True
+        return adaptor.get_X_array(None, varmask)
 
     def test_anndata_default(self):
         """Test an anndata adaptor with its default diffexp algorithm (diffexp_generic)"""
@@ -80,3 +88,64 @@ class DiffExpTest(unittest.TestCase):
         # run it directly
         results = diffexp_generic.diffexp_ttest(adaptor, maskA, maskB, 10)
         self.check_1_10_2_10(results)
+
+    def test_cxg_sparse(self):
+        self.sparse_diffexp(False)
+
+    def test_cxg_sparse_col_shift(self):
+        self.sparse_diffexp(True)
+
+    def sparse_diffexp(self, apply_col_shift):
+        with tempfile.TemporaryDirectory() as dirname:
+            # create a sparse matrix
+            h5adfile = os.path.join(dirname, "sparse.h5ad")
+            create_test_h5ad(h5adfile, 2000, 2000, 10, apply_col_shift)
+            adaptor_anndata = self.load_dataset(h5adfile, extra=dict(embeddings__names=[]))
+            adata = adaptor_anndata.data
+
+            sparsename = os.path.join(dirname, "sparse.cxg")
+            write_cxg(adata=adata, container=sparsename, title="sparse", sparse_threshold=11)
+            adaptor_sparse = self.load_dataset(sparsename)
+            assert adaptor_sparse.open_array("X").schema.sparse
+            assert adaptor_sparse.has_array("X_col_shift") == apply_col_shift
+
+            densename = os.path.join(dirname, "dense.cxg")
+            write_cxg(adata=adata, container=densename, title="dense", sparse_threshold=0)
+            adaptor_dense = self.load_dataset(densename)
+            assert not adaptor_dense.open_array("X").schema.sparse
+            assert not adaptor_dense.has_array("X_col_shift")
+
+            maskA = self.get_mask(adaptor_anndata, 1, 10)
+            maskB = self.get_mask(adaptor_anndata, 2, 10)
+
+            diffexp_results_anndata = diffexp_generic.diffexp_ttest(adaptor_anndata, maskA, maskB, 10)
+            diffexp_results_sparse = diffexp_cxg.diffexp_ttest(adaptor_sparse, maskA, maskB, 10)
+            diffexp_results_dense = diffexp_cxg.diffexp_ttest(adaptor_dense, maskA, maskB, 10)
+
+            self.compare_diffexp_results(diffexp_results_anndata, diffexp_results_sparse)
+            self.compare_diffexp_results(diffexp_results_anndata, diffexp_results_dense)
+
+            topcols = np.array([x[0] for x in diffexp_results_anndata])
+            cols_anndata = self.get_X_col(adaptor_anndata, topcols)
+            cols_sparse = self.get_X_col(adaptor_sparse, topcols)
+            cols_dense = self.get_X_col(adaptor_dense, topcols)
+            assert cols_anndata.shape[0] == adaptor_sparse.get_shape()[0]
+            assert cols_anndata.shape[1] == len(diffexp_results_anndata)
+
+            def convert(mat, cols):
+                return decode_matrix_fbs(encode_matrix_fbs(mat, col_idx=cols)).to_numpy()
+
+            cols_anndata = convert(cols_anndata, topcols)
+            cols_sparse = convert(cols_sparse, topcols)
+            cols_dense = convert(cols_dense, topcols)
+
+            x = adaptor_sparse.get_X_array()
+            assert x.shape == adaptor_sparse.get_shape()
+
+            for row in range(cols_anndata.shape[0]):
+                for col in range(cols_anndata.shape[1]):
+                    vanndata = cols_anndata[row][col]
+                    vsparse = cols_sparse[row][col]
+                    vdense = cols_dense[row][col]
+                    self.assertTrue(np.isclose(vanndata, vsparse, 1e-6, 1e-6))
+                    self.assertTrue(np.isclose(vanndata, vdense, 1e-6, 1e-6))