sparse column shift encoding. (#1502)

Many of our matrices are log normalized, which tends to eliminate
the number of non zero values (if there were any).  This prevents
the matrix from being stored as a sparse matrix.  The solution here
is to use a simple transformation to make it sparse again.  The most
common value from each column is subtracted from that column.  These
values that were subtracted are saved in an array called X_col_shift.

The cellxgene code needs to understand how to undo the transformation when
operating over the X matrix.

- added script to create a synthetic dataset for testing
- added a script to convert an existing CXG dataset to a sparse CXG dataset

server/compute/diffexp_cxg.py

@@ -107,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),

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
@@ -205,7 +209,7 @@ def write_cxg(
     log(1, "\t...embeddings created")
 
     # X matrix
-    save_X(container, adata, ctx, sparse_threshold)
+    save_X(container, adata.X, ctx, sparse_threshold)
     log(1, "\t...X created")
 
 
@@ -437,72 +441,146 @@ def create_X(X_name, shape, is_sparse):
         tiledb.DenseArray.create(X_name, schema)
 
 
-def check_sparse_X(adata, sparse_threshold):
-    shape = adata.X.shape
+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 = adata.X[row:lim, :]
+        a = xdata[row:lim, :]
         if type(a) is not np.ndarray:
             a = a.toarray()
         nnz += np.count_nonzero(a)
         if nnz > maxnnz:
-            return (nnz, lim * shape[1])
-        log(2, "\t...rows", row, "to", lim, "nnz", nnz, "nnz percent %5.2f%%" % (100 * nnz / (lim * shape[1])))
+            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])))
 
-    return (nnz, shape[0] * shape[1])
+    is_sparse = (100.0 * nnz / (shape[0] * shape[1])) < sparse_threshold
+    return (is_sparse, nnz, shape[0] * shape[1])
 
 
-def save_X(container, adata, ctx, sparse_threshold):
+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
 
+    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:
-        nnz, nelem = check_sparse_X(adata, sparse_threshold)
+        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...shape:", str(shape), "non-zeros percent (estimate): %6.2f" % percent)
+            log(1, "\t...sparse=", is_sparse, "non-zeros percent (estimate): %6.2f" % percent)
         else:
-            log(1, "\t...shape:", str(shape), "non-zeros:", nnz, "percent: %6.2f" % percent)
+            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)
     if is_sparse:
-        log(1, "\t...output X as sparse matrix")
+        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 = adata.X[row:lim, :]
+                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", row, "to", lim)
-            tiledb.consolidate(X_name, ctx=ctx)
+                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 = adata.X[row:lim, :]
+                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}")
@@ -251,10 +254,18 @@ class CxgAdaptor(DataAdaptor):
                 data = X[:, :]
             else:
                 data = X.multi_index[obs_items, var_items]
-            nrows, obsindices = self.__remap_indices(X.shape[0], obs_mask, data["obs"])
-            ncols, varindices = self.__remap_indices(X.shape[1], var_mask, data["var"])
+
+            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:

server/test/__init__.py

@@ -60,7 +60,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,6 +74,7 @@ def app_config(data_locator, backed=False):
     }
     config = AppConfig()
     config.update(**args)
+    config.update(**extra)
     config.complete_config()
     return config
 

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/test_diffexp.py

@@ -4,9 +4,10 @@ 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
 import tempfile
-import anndata
 import os
 
 
@@ -14,8 +15,8 @@ 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):
-        config = app_config(path)
+    def load_dataset(self, path, extra={}):
+        config = app_config(path, extra=extra)
         loader = MatrixDataLoader(path)
         adaptor = loader.open(config)
         return adaptor
@@ -28,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 = [
@@ -42,12 +51,7 @@ 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.assertTrue(np.isclose(result[1], expect[1], 1e-6, 1e-6))
-            self.assertTrue(np.isclose(result[2], expect[2], 1e-6, 1e-6))
-            self.assertTrue(np.isclose(result[3], expect[3], 1e-6, 1e-6))
+        self.compare_diffexp_results(results, expects)
 
     def get_X_col(self, adaptor, cols):
         varmask = np.zeros(adaptor.get_shape()[1], dtype=bool)
@@ -86,33 +90,62 @@ class DiffExpTest(unittest.TestCase):
         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")
-            densename = os.path.join(dirname, "dense.cxg")
-            source_h5ad = anndata.read_h5ad(f"{PROJECT_ROOT}/example-dataset/pbmc3k.h5ad")
-            # create a cxg sparse array
-            write_cxg(adata=source_h5ad, container=sparsename, title="pbmc3k", sparse_threshold=100)
-            write_cxg(adata=source_h5ad, container=densename, title="pbmc3k", sparse_threshold=0)
+            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")
 
-            col_results = []
-            for adaptor in (adaptor_sparse, adaptor_dense):
-                maskA = self.get_mask(adaptor, 1, 10)
-                maskB = self.get_mask(adaptor, 2, 10)
-                diffexp_results = diffexp_cxg.diffexp_ttest(adaptor, maskA, maskB, 10)
-                self.check_1_10_2_10(diffexp_results)
-                topcols = [x[0] for x in diffexp_results]
-                cols = self.get_X_col(adaptor, topcols)
-                assert cols.shape[0] == adaptor.get_shape()[0]
-                assert cols.shape[1] == len(diffexp_results)
-                col_results.append(cols)
+            maskA = self.get_mask(adaptor_anndata, 1, 10)
+            maskB = self.get_mask(adaptor_anndata, 2, 10)
 
-                x = adaptor.get_X_array()
-                print(x)
+            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)
 
-        for row in range(col_results[0].shape[0]):
-            for col in range(col_results[0].shape[1]):
-                sval = col_results[0][row][col]
-                dval = col_results[1][row][col]
-                self.assertTrue(np.isclose(sval, dval, 1e-6, 1e-6))
+            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))