From 994c20c0945ed0ede098e67df5465e7512dce584 Mon Sep 17 00:00:00 2001 From: maniarathi Date: Mon, 17 Aug 2020 17:28:29 -0700 Subject: [PATCH] Move cxgtool into CLI and modularize conversion functions (#1701) --- server/cli/convert_to_cxg.py | 131 ++++ server/common/annotations/hosted_tiledb.py | 17 +- server/common/app_config.py | 2 +- server/common/corpora.py | 20 +- server/common/utils/corpora_constants.py | 16 + server/common/utils/cxg_constants.py | 4 + server/common/utils/cxg_generation_utils.py | 179 +++++ server/common/utils/type_conversion_utils.py | 11 + server/converters/cxgtool.py | 669 ------------------ server/converters/h5ad_data_file.py | 250 +++++++ server/test/unit/common/test_api.py | 6 +- server/test/unit/common/test_corpora.py | 9 +- .../common/utils/test_cxg_generation_utils.py | 149 ++++ .../utils/test_type_conversion_utils.py | 18 +- server/test/unit/compute/test_diffexp_cxg.py | 33 +- server/test/unit/converters/test_cxgtool.py | 41 -- .../unit/converters/test_h5ad_data_file.py | 235 ++++++ 17 files changed, 1032 insertions(+), 758 deletions(-) create mode 100644 server/cli/convert_to_cxg.py create mode 100644 server/common/utils/corpora_constants.py create mode 100644 server/common/utils/cxg_constants.py create mode 100644 server/common/utils/cxg_generation_utils.py delete mode 100644 server/converters/cxgtool.py create mode 100644 server/converters/h5ad_data_file.py create mode 100644 server/test/unit/common/utils/test_cxg_generation_utils.py delete mode 100644 server/test/unit/converters/test_cxgtool.py create mode 100644 server/test/unit/converters/test_h5ad_data_file.py diff --git a/server/cli/convert_to_cxg.py b/server/cli/convert_to_cxg.py new file mode 100644 index 00000000..e8d37936 --- /dev/null +++ b/server/cli/convert_to_cxg.py @@ -0,0 +1,131 @@ +from os import path + +import click + +from server.converters.h5ad_data_file import H5ADDataFile + + +@click.command( + name="convert", + short_help="Converts an H5AD dataset to the CXG format.", + help="Converts an H5AD dataset to the CXG format. The CXG format is a cellxgene-private data format " + "that has performance and access characteristics amenable to a multi-dataset, multi-user serving " + "environment. You will be able to launch the cellxgene using the `cellxgene launch` command as " + "usually with the generated CXG file.", +) +@click.argument( + "input-file", + nargs=1, + help="Path to the H5AD input file to be converted.", + type=click.Path(exists=True, dir_okay=False), +) +@click.option( + "-o", + "--output-dir", + help="Name of the output CXG directory. If not provided, will default to be the input filename with a " + "CXG extension.", +) +@click.option( + "-b", + "--backed", + help="When true, loads the H5AD in file backed mode. This will cause the conversion to be slower, " + "but will use less memory.", + default=False, + show_default=True, + is_flag=True, +) +@click.option( + "-t", + "--title", + help="Human readable dataset title that will be included as metadata about the CXG file. If omitted, " + "the dataset title will be the filename.", +) +@click.option( + "-a", + "--about", + help="A fully qualified URL that provides more information about the dataset and will be included as " + "metadata about the CXG file.", +) +@click.option( + "-s", + "--sparse-threshold", + help="If the dataset's percent of non-zero values falls belows the specified threshold, then the X " + "array of the dataset will be sparse. Since the default value is 0.0, the default will be to " + "convert to dense array.", + default=0.0, + show_default=True, +) +@click.option("--obs-names", + help="Name to a column in the obs dataframe that will be used as the index for the dataframe instead of " + "the one designated by the dataframe generated-index.") +@click.option("--var-names", + help="Name to a column in the var dataframe that will be used as the index for the dataframe instead of " + "the one designated by the dataframe generated-index.") +@click.option( + "--disable-custom-colors", + help="When set, conversion process will not extract scanpy-compatible category colors from the H5AD file.", + default=False, + show_default=True, + is_flag=True, +) +@click.option( + "--disable-corpora-schema", + "When set, conversion process will neither extract nor store Corpora schema information. See " + "https://github.com/chanzuckerberg/corpora-data-portal/blob/main/backend/schema/corpora_schema.md for " + "more information.", + default=False, + show_default=True, + is_flag=True, +) +@click.option( + "--overwrite", + help="When set to true, will overwrite the output file if the output file already exists.", + default=False, + show_default=True, + is_flag=True, +) +@click.option("-v", "--verbose", count=True) +@click.help_option("--help", "-h", help="Show this message and exit.") +def convert_to_cxg( + input_file, + output_directory, + backed, + title, + about, + sparse_threshold, + obs_names, + var_names, + disable_custom_colors, + disable_corpora_schema, + should_overwrite, +): + """ + Convert a dataset file into CXG. + """ + + h5ad_data_file = H5ADDataFile(input_file, backed, title, about, obs_names, var_names, + use_corpora_schema=not disable_corpora_schema) + + # Get the directory that will hold all the CXG files + cxg_output_container = get_output_directory(input_file, output_directory, should_overwrite) + + h5ad_data_file.to_cxg(cxg_output_container, sparse_threshold, + convert_anndata_colors_to_cxg_colors=not disable_custom_colors) + + +def get_output_directory(input_filename, output_directory, should_overwrite): + """ + Get the name of the CXG output directory to be created/populated during the dataset conversion. + """ + + if not path.isdir(output_directory) or (path.isdir(output_directory) and should_overwrite): + if output_directory.endswith(".cxg"): + return output_directory + return output_directory + ".cxg" + if path.isdir(output_directory) and not should_overwrite: + raise click.BadParameter( + f"Output directory {output_directory} already exists. If you'd like to overwrite, then run the command " + f"with the --overwrite flag." + ) + + return path.splitext(input_filename)[1] + ".cxg" diff --git a/server/common/annotations/hosted_tiledb.py b/server/common/annotations/hosted_tiledb.py index 05b017a4..16ac7b29 100644 --- a/server/common/annotations/hosted_tiledb.py +++ b/server/common/annotations/hosted_tiledb.py @@ -8,7 +8,9 @@ import tiledb from flask import current_app from server.common.annotations.annotations import Annotations -from server.converters.cxgtool import sanitize_keys, generate_schema_hints_and_convert_value_types, cxg_dtype +from server.common.errors import AnnotationCategoryNameError +from server.common.utils.sanitization_utils import sanitize_values_in_list +from server.common.utils.type_conversion_utils import get_dtypes_and_schemas_of_dataframe, get_dtype_of_array from server.db.cellxgene_orm import CellxGeneDataset, Annotation @@ -21,7 +23,12 @@ class AnnotationsHostedTileDB(Annotations): self.directory_path = directory_path def check_category_names(self, df): - sanitize_keys(df.keys().to_list(), False) + original_category_names = df.keys().to_list() + sanitized_category_names = set(sanitize_values_in_list(original_category_names).values()) + unsanitary_original_category_names = set(original_category_names).difference(sanitized_category_names) + if unsanitary_original_category_names: + raise AnnotationCategoryNameError( + f"{unsanitary_original_category_names} are not valid category names, please resubmit") def is_safe_collection_name(self, name): """ @@ -94,19 +101,19 @@ class AnnotationsHostedTileDB(Annotations): pass else: os.makedirs(uri, exist_ok=True) - schema_hints, values = generate_schema_hints_and_convert_value_types(df) + _, dataframe_schema_type_hints = get_dtypes_and_schemas_of_dataframe(df) annotation = Annotation( tiledb_uri=uri, user_id=user_id, dataset_id=str(dataset_id), - schema_hints=json.dumps(schema_hints) + schema_hints=json.dumps(dataframe_schema_type_hints) ) if not df.empty: self.check_category_names(df) # convert to tiledb datatypes for col in df: - df[col] = df[col].astype(cxg_dtype(df[col])) + df[col] = df[col].astype(get_dtype_of_array(df[col])) tiledb.from_pandas(uri, df) self.db.session.add(annotation) diff --git a/server/common/app_config.py b/server/common/app_config.py index e262bc49..df7a0f96 100644 --- a/server/common/app_config.py +++ b/server/common/app_config.py @@ -768,7 +768,7 @@ class DatasetConfig(BaseConfig): self.user_annotations__ontology__obo_location = dc["user_annotations"]["ontology"]["obo_location"] self.user_annotations__hosted_tiledb_array__db_uri = dc["user_annotations"]["hosted_tiledb_array"]["db_uri"] self.user_annotations__hosted_tiledb_array__hosted_file_directory = \ - dc["user_annotations"]["hosted_tiledb_array"]["hosted_file_directory"] # noqa E501 + dc["user_annotations"]["hosted_tiledb_array"]["hosted_file_directory"] # noqa E501 self.embeddings__names = dc["embeddings"]["names"] self.embeddings__enable_reembedding = dc["embeddings"]["enable_reembedding"] diff --git a/server/common/corpora.py b/server/common/corpora.py index 9f48873b..d616f75b 100644 --- a/server/common/corpora.py +++ b/server/common/corpora.py @@ -9,6 +9,7 @@ import collections import json from server.cli.upgrade import validate_version_str +from server.common.utils.corpora_constants import CorporaConstants def corpora_get_versions_from_anndata(adata): @@ -56,26 +57,13 @@ def corpora_get_props_from_anndata(adata): if not version_is_supported: raise ValueError("Unsupported Corpora schema version") - required_simple_fields = [ - "version", - "title", - "layer_descriptions", - "organism", - "organism_ontology_term_id", - "project_name", - "project_description", - ] - # Spec says some values encoded as JSON due to the inability of AnnData to store complex types. - required_json_fields = ["contributors", "project_links"] - optional_simple_fields = ["preprint_doi", "publication_doi", "default_embedding", "default_field", "tags"] - corpora_props = {} - for key in required_simple_fields: + for key in CorporaConstants.REQUIRED_SIMPLE_METADATA_FIELDS: if key not in adata.uns: raise KeyError(f"missing Corpora schema field {key}") corpora_props[key] = adata.uns[key] - for key in required_json_fields: + for key in CorporaConstants.REQUIRED_JSON_ENCODED_METADATA_FIELD: if key not in adata.uns: raise KeyError(f"missing Corpora schema field {key}") try: @@ -83,7 +71,7 @@ def corpora_get_props_from_anndata(adata): except json.JSONDecodeError: raise json.JSONDecodeError(f"Corpora schema field {key} is expected to be a valid JSON string") - for key in optional_simple_fields: + for key in CorporaConstants.OPTIONAL_SIMPLE_METADATA_FIELDS: if key in adata.uns: corpora_props[key] = adata.uns[key] diff --git a/server/common/utils/corpora_constants.py b/server/common/utils/corpora_constants.py new file mode 100644 index 00000000..c0c98168 --- /dev/null +++ b/server/common/utils/corpora_constants.py @@ -0,0 +1,16 @@ +class CorporaConstants(object): + REQUIRED_SIMPLE_METADATA_FIELDS = [ + "version", + "title", + "layer_descriptions", + "organism", + "organism_ontology_term_id", + "project_name", + "project_description", + ] + + # The Corpora specification requires some values encoded as JSON due to the inability of AnnData to store complex + # types. + REQUIRED_JSON_ENCODED_METADATA_FIELD = ["contributors", "project_links"] + + OPTIONAL_SIMPLE_METADATA_FIELDS = ["preprint_doi", "publication_doi", "default_embedding", "default_field", "tags"] diff --git a/server/common/utils/cxg_constants.py b/server/common/utils/cxg_constants.py new file mode 100644 index 00000000..5982ec15 --- /dev/null +++ b/server/common/utils/cxg_constants.py @@ -0,0 +1,4 @@ +class CxgConstants(object): + # The CXG container version number. Must be a semver string (major.minor.patch) + # DO NOT UPDATE THIS WITHOUT ALSO UPDATING CXG SPECIFICATION. + CXG_VERSION = "0.2.0" diff --git a/server/common/utils/cxg_generation_utils.py b/server/common/utils/cxg_generation_utils.py new file mode 100644 index 00000000..f336c0ae --- /dev/null +++ b/server/common/utils/cxg_generation_utils.py @@ -0,0 +1,179 @@ +import json + +import numpy as np +import tiledb + +from server.common.utils.type_conversion_utils import get_dtype_of_array, get_dtype_and_schema_of_array + + +def convert_dictionary_to_cxg_group(cxg_container, metadata_dict, group_metadata_name="cxg_group_metadata"): + """ + Saves the contents of the dictionary to the CXG output directory specified. + + This function is primarily used to save metadata about a dataset to the CXG directory. At some point, tiledb will + have support for metadata on groups at which point the utility of this function should be revisited. Until such + feature exists, this function create an empty array and annotate that array. + + For more information, visit https://github.com/TileDB-Inc/TileDB-Py/issues/254. + """ + + array_name = f"{cxg_container}/{group_metadata_name}" + + # Because TileDB does not allow one to attach metadata directly to a CXG group, we need to have a workaround + # where we create an empty array and attached the metadata onto to this empty array. Below we construct this empty + # array. + tiledb.from_numpy(array_name, np.zeros((1,))) + + with tiledb.DenseArray(array_name, mode="w") as metadata_array: + for key, value in metadata_dict.items(): + metadata_array.meta[key] = value + + +def convert_dataframe_to_cxg_array(cxg_container, dataframe_name, dataframe, index_column_name, ctx): + """ + Saves the contents of the dataframe to the CXG output directory specified. + + Current access patterns are oriented toward reading very large slices of the dataframe, one attribute at a time. + Attribute data also tends to be (often) repetitive (bools, categories, strings). Given this, we use a large tile + size (1000) and very aggressive compression levels. + """ + + def create_dataframe_array(array_name, dataframe): + tiledb_filter = tiledb.FilterList( + [ + # Attempt aggressive compression as many of these dataframes are very repetitive strings, bools and + # other non-float data. + tiledb.ZstdFilter(level=22), + ] + ) + attrs = [ + tiledb.Attr(name=column, dtype=get_dtype_of_array(dataframe[column]), filters=tiledb_filter) + for column in dataframe + ] + domain = tiledb.Domain( + tiledb.Dim(domain=(0, dataframe.shape[0] - 1), tile=min(dataframe.shape[0], 1000), dtype=np.uint32) + ) + schema = tiledb.ArraySchema( + domain=domain, sparse=False, attrs=attrs, cell_order="row-major", tile_order="row-major" + ) + tiledb.DenseArray.create(array_name, schema) + + array_name = f"{cxg_container}/{dataframe_name}" + + create_dataframe_array(array_name, dataframe) + + with tiledb.DenseArray(array_name, mode="w", ctx=ctx) as array: + value = {} + schema_hints = {} + for column_name, column_values in dataframe.items(): + dtype, hints = get_dtype_and_schema_of_array(column_values) + + value[column_name] = column_values.to_numpy(dtype=dtype) + if hints: + schema_hints.update({column_name: hints}) + + schema_hints.update({"index": index_column_name}) + array[:] = value + array.meta["cxg_schema"] = json.dumps(schema_hints) + + tiledb.consolidate(array_name, ctx=ctx) + + +def convert_ndarray_to_cxg_dense_array(ndarray_name, ndarray, ctx): + """ + Saves contents of ndarray to the CXG output directory specified. + + Generally this function is used to convert dataset embeddings. Because embeddings are typically accessed with + very large slices (or all of the embedding), they do not benefit from overly aggressive compression due to their + format. Given this, we use a large tile size (1000) but only default compression level. + """ + + def create_ndarray_array(ndarray_name, ndarray): + filters = tiledb.FilterList([tiledb.ZstdFilter()]) + attrs = [tiledb.Attr(dtype=ndarray.dtype, filters=filters)] + dimensions = [ + tiledb.Dim( + domain=(0, ndarray.shape[dimension] - 1), tile=min(ndarray.shape[dimension], 1000), dtype=np.uint32 + ) + for dimension in range(ndarray.ndim) + ] + domain = tiledb.Domain(*dimensions) + schema = tiledb.ArraySchema( + domain=domain, sparse=False, attrs=attrs, capacity=1_000_000, cell_order="row-major", tile_order="row-major" + ) + tiledb.DenseArray.create(ndarray_name, schema) + + create_ndarray_array(ndarray_name, ndarray) + + with tiledb.DenseArray(ndarray_name, mode="w", ctx=ctx) as array: + array[:] = ndarray + + tiledb.consolidate(ndarray_name, ctx=ctx) + + +def convert_matrix_to_cxg_array( + matrix_name, matrix, encode_as_sparse_array, ctx, column_shift_for_sparse_encoding=None +): + """ + Converts a numpy array matrix into a TileDB SparseArray of DenseArray based on whether `encode_as_sparse_array` + is true or not. Note that when the matrix is encoded as a SparseArray, it only writes the values that are + nonzero. This means that if you count the number of elements in the SparseArray, it will not equal the total + number of elements in the matrix, only the number of nonzero elements. + + Furthermore, if the `column_shift_for_sparse_encoding` matrix is not None, this function will subtract the sparse + encoding from the original given matrix and as previously stated, only write the nonzero values to the TileDB + SparseArray. + """ + + def create_matrix_array(matrix_name, number_of_rows, number_of_columns, encode_as_sparse_array): + filters = tiledb.FilterList([tiledb.ZstdFilter()]) + attrs = [tiledb.Attr(dtype=np.float32, filters=filters)] + if encode_as_sparse_array: + domain = tiledb.Domain( + tiledb.Dim(name="obs", domain=(0, number_of_rows - 1), tile=min(number_of_rows, 512), dtype=np.uint32), + tiledb.Dim( + name="var", domain=(0, number_of_columns - 1), tile=min(number_of_columns, 2048), dtype=np.uint32 + ), + ) + else: + domain = tiledb.Domain( + tiledb.Dim(name="obs", domain=(0, number_of_rows - 1), tile=min(number_of_rows, 50), dtype=np.uint32), + tiledb.Dim( + name="var", domain=(0, number_of_columns - 1), tile=min(number_of_columns, 100), dtype=np.uint32 + ), + ) + schema = tiledb.ArraySchema( + domain=domain, sparse=encode_as_sparse_array, attrs=attrs, cell_order="row-major", tile_order="col-major" + ) + if encode_as_sparse_array: + tiledb.SparseArray.create(matrix_name, schema) + else: + tiledb.DenseArray.create(matrix_name, schema) + + number_of_rows = matrix.shape[0] + number_of_columns = matrix.shape[1] + stride = min(int(np.power(10, np.around(np.log10(1e9 / number_of_columns)))), 10_000) + + create_matrix_array(matrix_name, number_of_rows, number_of_columns, encode_as_sparse_array) + + if encode_as_sparse_array: + with tiledb.SparseArray(matrix_name, mode="w", ctx=ctx) as array: + for start_row_index in range(0, number_of_rows, stride): + end_row_index = min(start_row_index + stride, number_of_rows) + matrix_subset = matrix[start_row_index:end_row_index, :] + if not isinstance(matrix_subset, np.ndarray): + matrix_subset = matrix_subset.toarray() + if column_shift_for_sparse_encoding is not None: + matrix_subset = matrix_subset - column_shift_for_sparse_encoding + indices = np.nonzero(matrix_subset) + trow = indices[0] + start_row_index + array[trow, indices[1]] = matrix_subset[indices[0], indices[1]] + + else: + with tiledb.DenseArray(matrix_name, mode="w", ctx=ctx) as array: + for start_row_index in range(0, number_of_rows, stride): + end_row_index = min(start_row_index + stride, number_of_rows) + matrix_subset = matrix[start_row_index:end_row_index, :] + if not isinstance(matrix_subset, np.ndarray): + matrix_subset = matrix_subset.toarray() + array[start_row_index:end_row_index, :] = matrix_subset diff --git a/server/common/utils/type_conversion_utils.py b/server/common/utils/type_conversion_utils.py index 8b467b47..eca8b938 100644 --- a/server/common/utils/type_conversion_utils.py +++ b/server/common/utils/type_conversion_utils.py @@ -4,6 +4,17 @@ import numpy as np import pandas as pd +def get_dtypes_and_schemas_of_dataframe(dataframe: pd.DataFrame): + dtypes_by_column_name = {} + schema_type_hints_by_column_name = {} + + for column_name, column_values in dataframe.items(): + dtypes_by_column_name[column_name], schema_type_hints_by_column_name[column_name] = \ + get_dtype_and_schema_of_array(column_values) + + return dtypes_by_column_name, schema_type_hints_by_column_name + + def get_dtype_of_array(array: pd.Series): return get_dtype_and_schema_of_array(array)[0] diff --git a/server/converters/cxgtool.py b/server/converters/cxgtool.py deleted file mode 100644 index 5116c4e6..00000000 --- a/server/converters/cxgtool.py +++ /dev/null @@ -1,669 +0,0 @@ -""" -This program converts an [AnnData H5AD](https://anndata.readthedocs.io/en/stable/) -into a cellxgene TileDB structure, aka a [CXG](../../dev_docs/cxg.md). - -IF YOU UPDATE THIS FILE, IN ANY WAY THAT MODIFIES THE CXG FORMAT or CONTENTS, -YOU MUST UPDATE THE CXG SPECIFICATION and VERSION NUMBER. -""" -import re -import anndata -import tiledb -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, AnnotationCategoryNameError -from server.common.corpora import ( - corpora_get_props_from_anndata, - corpora_get_versions_from_anndata, - corpora_is_version_supported, -) - - -# the CXG container version number. Must be a semver string (major.minor.patch) -# DO NOT UPDATE THIS WITHOUT ALSO UPDATING THE CXG SPECIFICATION. -CXG_VERSION = "0.2.0" - -# log_level must have a default -log_level = 3 - - -def log(level, *args): - global log_level - if log_level and level <= log_level: - print(*args) - - -def main(): - parser = argparse.ArgumentParser() - parser.add_argument("h5ad", nargs="?", help="H5AD file name") - parser.add_argument( - "--backed", action="store_true", help="loaded in file backed mode. Will be slower, but use less memory." - ) - parser.add_argument( - "--disable-custom-colors", - action="store_true", - default=False, - help="Do not extract scanpy-compatible category colors from h5ad file.", - ) - parser.add_argument( - "--obs-names", help="Name of annotation to use for observations. If not specified, will use the obs index." - ) - parser.add_argument( - "--var-names", help="Name of annotation to use for variables. If not specified, will use the var index." - ) - parser.add_argument("--verbose", "-v", action="count", default=0, help="verbose output") - parser.add_argument("--title", help="Human readable dataset title. If omitted, will use filename") - parser.add_argument( - "--about", - metavar="", - 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", - ) - parser.add_argument( - "--disable-corpora", - action="store_true", - default=False, - help="Disable extraction and storing of Corpora schema information.", - ) - args = parser.parse_args() - - global log_level - log_level = args.verbose - - adata = anndata.read_h5ad(args.h5ad, backed="r" if args.backed else None) - log(1, f"{basename(args.h5ad)} loaded...") - - basefname = splitext(basename(args.h5ad))[0] - out = args.out if args.out is not None else basefname - container = out if splitext(out)[1] == ".cxg" else out + ".cxg" - - corpora_props = load_corpora_props(args, adata) if not args.disable_corpora else None - cxg_group_metadata = create_cxg_group_metadata( - adata, - basefname, - title=args.title, - about=args.about, - corpora_props=corpora_props, - extract_colors=not args.disable_custom_colors, - ) - - write_cxg( - adata, - container, - cxg_group_metadata=cxg_group_metadata, - var_names=args.var_names, - obs_names=args.obs_names, - sparse_threshold=args.sparse_threshold, - ) - - log(1, "done") - - -def write_cxg(adata, container, cxg_group_metadata, var_names=None, obs_names=None, 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: - raise ValueError("Observation index is not unique - unable to convert.") - - """ - TileDB bug TileDB-Inc/TileDB#1575 requires that we sanitize all column names - prior to saving. This can be reverted when the bug is fixed. - """ - log(0, "Warning: sanitizing all dataframe column names.") - clean_all_column_names(adata) - - ctx = tiledb.Ctx( - { - "sm.num_reader_threads": 32, - "sm.num_writer_threads": 32, - "sm.consolidation.buffer_size": 1 * 1024 * 1024 * 1024, - } - ) - - tiledb.group_create(container, ctx=ctx) - log(1, f"\t...group created, with name {container}") - - # dataset metadata - save_metadata(container, cxg_group_metadata) - log(1, "\t...dataset metadata saved") - - # var/gene dataframe - save_dataframe(container, "var", adata.var, var_names, ctx=ctx) - log(1, "\t...var dataframe created") - - # obs/cell dataframe - save_dataframe(container, "obs", adata.obs, obs_names, ctx=ctx) - log(1, "\t...obs dataframe created") - - # embeddings - e_container = f"{container}/emb" - tiledb.group_create(e_container, ctx=ctx) - save_embeddings(e_container, adata, ctx) - log(1, "\t...embeddings created") - - # X matrix - save_X(container, adata.X, ctx, sparse_threshold) - log(1, "\t...X created") - - -""" -TODO: the code used to handle type inferencing should not be duplicated between -this tool and the server/common/utils code. When this tool is merged into -the cellxgene CLI, consolidate. -""" - - -def dtype_to_schema(dtype): - if dtype == np.float32: - return (np.float32, {}) - elif dtype == np.int32: - return (np.int32, {}) - elif dtype == np.bool_: - return (np.uint8, {"type": "boolean"}) - elif dtype == np.str: - return (np.unicode, {"type": "string"}) - elif dtype == "category": - typ, hint = cxg_type(dtype.categories) - return (typ, {"type": "categorical", "categories": dtype.categories.tolist()}) - else: - raise TypeError(f"Annotations of type {dtype} are unsupported.") - - -def _can_cast_to_float32(array): - if array.dtype.kind == "f": - # force downcast for all floats - return True - return False - - -def _can_cast_to_int32(array): - if array.dtype.kind in ["i", "u"]: - if np.can_cast(array.dtype, np.int32): - return True - ii32 = np.iinfo(np.int32) - if array.min() >= ii32.min and array.max() <= ii32.max: - return True - return False - - -def cxg_type(array): - try: - return dtype_to_schema(array.dtype) - except TypeError: - dtype = array.dtype - data_kind = dtype.kind - if _can_cast_to_float32(array): - return (np.float32, {}) - elif _can_cast_to_int32(array): - return (np.int32, {}) - elif data_kind == "O" and dtype == "object": - return (np.unicode, {"type": "string"}) - else: - raise TypeError(f"Annotations of type {dtype} are unsupported.") - - -def cxg_dtype(array): - return cxg_type(array)[0] - - -def create_dataframe(name, df, ctx): - """ - Current access patterns are oriented toward reading very large slices of - the dataframe, one attribute at a time. Attribute data also tends to be - (often) repetitive (bools, categories, strings). - Given this, we use: - * a large tile size (1000) - * very aggressive compression levels - """ - filter = tiledb.FilterList( - [ - # attempt aggressive compression as many of these dataframes are very repetitive - # strings, bools and other non-float data. - tiledb.ZstdFilter(level=22), - ] - ) - attrs = [tiledb.Attr(name=col, dtype=cxg_dtype(df[col]), filters=filter) for col in df] - domain = tiledb.Domain(tiledb.Dim(domain=(0, df.shape[0] - 1), tile=min(df.shape[0], 1000), dtype=np.uint32)) - schema = tiledb.ArraySchema( - domain=domain, sparse=False, attrs=attrs, cell_order="row-major", tile_order="row-major" - ) - tiledb.DenseArray.create(name, schema) - - -def create_unique_column_name(df_cols, col_name_prefix): - """ - given the columns of a dataframe, and a name prefix, return a column name which - does not exist in the dataframe, AND which is prefixed by `prefix` - - The approach is to append a numeric suffix, starting at zero and increasing by - one, until an unused name is found (eg, prefix_0, prefix_1, ...). - """ - suffix = 0 - while f"{col_name_prefix}{suffix}" in df_cols: - suffix += 1 - return f"{col_name_prefix}{suffix}" - - -def alias_index_col(df, df_name, index_col_name): - """ - We rely in the existance of a unique, human-readable index for - any dataframe (eg, var is typically gene name, obs the cell name). - The user can specify these via the --obs-names and --var-names config. - If they are not specified, use the existing index to create them, giving - the resulting column a unique name (eg, "name"). - - In both cases, enforce that the result is unique, and communicate the - index column name via the 'index' field in the schema hints. - """ - if index_col_name is None: - if not df.index.is_unique: - raise KeyError( - f"Values in {df_name}.index must be unique. " - "Please prepare data to contain unique index values, or specify an " - "alternative with --{ax_name}-name." - ) - index_col_name = create_unique_column_name(df.columns, "name_") - # turn the index into a normal column - df.rename_axis(index_col_name, inplace=True) - df.reset_index(inplace=True) - - elif index_col_name in df.columns: - # User has specified alternative column for unique names, and it exists - if not df[index_col_name].is_unique: - raise KeyError( - f"Values in {df_name}.{index_col_name} must be unique. Please prepare data to contain unique values." - ) - - else: - raise KeyError(f"Annotation {index_col_name}, specified in --{df_name}-name, does not exist.") - - return (df, index_col_name) - - -def generate_schema_hints_and_convert_value_types(df): - value = {} - schema_hints = {} - for k, v in df.items(): - dtype, hints = cxg_type(v) - value[k] = v.to_numpy(dtype=dtype) - if hints: - schema_hints.update({k: hints}) - return schema_hints, value - - -def save_dataframe(container, name, df, index_col_name, ctx): - A_name = f"{container}/{name}" - (df, index_col_name) = alias_index_col(df, name, index_col_name) - create_dataframe(A_name, df, ctx=ctx) - with tiledb.DenseArray(A_name, mode="w", ctx=ctx) as A: - schema_hints, value = generate_schema_hints_and_convert_value_types(df) - schema_hints.update({"index": index_col_name}) - # convert all values in all cols to a numpy version of cxg datatypes, - # then store the contents in the tiledb array A - A[:] = value - A.meta["cxg_schema"] = json.dumps(schema_hints) - - tiledb.consolidate(A_name, ctx=ctx) - - -def create_emb(e_name, emb): - """ - Embeddings are typically accessed with very large slices (or all of the embedding), - and do not benefit from overly aggressive compression due to their format. Given - this, we use: - * large tile size (1000) - * default compression level - """ - filters = tiledb.FilterList([tiledb.ZstdFilter()]) - attrs = [tiledb.Attr(dtype=emb.dtype, filters=filters)] - 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)) - 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" - ) - tiledb.DenseArray.create(e_name, schema) - - -def is_valid_embedding(adata, name, arr): - """ return True if this layout data is a valid array for front-end presentation: - * ndarray, with shape (n_obs, >= 2), dtype float/int/uint - * follows ScanPy embedding naming conventions - * with all values finite or NaN (no +Inf or -Inf) - """ - is_valid = type(name) == str and name.startswith("X_") and len(name) > 2 - is_valid = is_valid and type(arr) == np.ndarray and arr.dtype.kind in "fiu" - is_valid = is_valid and arr.shape[0] == adata.n_obs and arr.shape[1] >= 2 - is_valid = is_valid and not np.any(np.isinf(arr)) and not np.all(np.isnan(arr)) - return is_valid - - -def save_embeddings(container, adata, ctx): - for (name, value) in adata.obsm.items(): - if is_valid_embedding(adata, name, value): - e_name = f"{container}/{name[2:]}" - create_emb(e_name, value) - with tiledb.DenseArray(e_name, mode="w", ctx=ctx) as A: - A[:] = value - tiledb.consolidate(e_name, ctx=ctx) - log(1, f"\t\t...{name} embedding created") - - -def create_X(X_name, shape, is_sparse): - """ - 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) 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)] - 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=is_sparse, attrs=attrs, cell_order="row-major", tile_order="col-major" - ) - if is_sparse: - tiledb.SparseArray.create(X_name, schema) - else: - tiledb.DenseArray.create(X_name, schema) - - -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 = 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) - 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): - """ - Save all dataset-wide metadata. This includes: - * CXG version - * dataset metadata, such as title and about link. - - Longer term, tiledb will have support for metadata on groups. Until - such feature exists, create an empty array and annotate that array. - - https://github.com/TileDB-Inc/TileDB-Py/issues/254 - """ - a_name = f"{container}/cxg_group_metadata" - with tiledb.from_numpy(a_name, np.zeros((1,))) as A: - pass - with tiledb.DenseArray(a_name, mode="w") as A: - for k, v in metadata_dict.items(): - A.meta[k] = v - - -def load_corpora_props(args, adata): - versions = corpora_get_versions_from_anndata(adata) - if versions is None: - return None - - [corpora_schema_version, corpora_encoding_version] = versions - corpora_props = corpora_get_props_from_anndata(adata) - version_is_supported = corpora_is_version_supported(corpora_schema_version, corpora_encoding_version) - if not version_is_supported or not corpora_props: - log(0, "ERROR: Unknown source file schema version is unsupported") - raise ValueError("Unsupported Corpora schema version") - - log(1, "FYI, file appears to be encoded using Corpora schema standards...") - if args.title is not None or args.about is not None: - log(0, "Warning: explicit specification of --title or --about will override Corpora schema fields.") - - return corpora_props - - -def create_cxg_group_metadata(adata, basefname, title=None, about=None, corpora_props=None, extract_colors=True): - - if corpora_props is not None: - # clobber encoding version to be OUR version, not the source H5AD encoding - corpora_props["version"].update({"corpora_encoding_version": CXG_VERSION}) - corpora_project_links = corpora_props.get("project_links", []) - corpora_about_link = next( - (link for link in corpora_project_links if (link.get("link_type", None) == "SUMMARY")), {} - ) - else: - corpora_about_link = {} - - title = title or corpora_about_link.get("link_name", basefname) - about = about or corpora_about_link.get("link_url") - - cxg_group_metadata = {"cxg_version": CXG_VERSION, "cxg_properties": json.dumps({"title": title, "about": about})} - if corpora_props is not None: - cxg_group_metadata.update({"corpora": json.dumps(corpora_props)}) - - if extract_colors: - try: - cxg_group_metadata["cxg_category_colors"] = json.dumps( - convert_anndata_category_colors_to_cxg_category_colors(adata) - ) - except ColorFormatException: - log( - 0, - "Warning: failed to extract colors from h5ad file! " - "Fix the h5ad file or rerun with --disable-custom-colors. See help for details.", - ) - - return cxg_group_metadata - - -def sanitize_keys(keys, update_keys=True): - """ - We need names to be safe to use as attribute names in tiledb. See: - TileDB-Inc/TileDB#1575 - TileDB-Inc/TileDB-Py#294 - This can be entirely removed once they add proper escaping. - - Args: list of keys - Returns: dict of {old_key: new_key, ...} - - Returned new keys will be both safe and unique. - - Masking out [~/.] and anything outside the ASCII range. - """ - p = re.compile(r"[^ -\.0-\[\]-\}]") - clean_keys = {k: p.sub("_", k) for k in keys} - - used_keys = set() - clean_unique_keys = {} - for k, v in clean_keys.items(): - if v not in used_keys: - used_keys.add(v) - clean_unique_keys[k] = v - continue - - # else, needs deduping. - counter = 1 - while True: - candidate_name = v + "-" + str(counter) - if candidate_name not in used_keys: - used_keys.add(candidate_name) - clean_unique_keys[k] = candidate_name - break - counter += 1 - - for k, v, in clean_unique_keys.items(): - if k != v: - if update_keys is False: - raise AnnotationCategoryNameError(f"{k} not a valid category name, please resubmit") - log(1, f"Renaming {k} to {v}") - return clean_unique_keys - - -def sanitize_df(df): - df.rename(columns=sanitize_keys(df.keys().tolist()), inplace=True) - - -def sanitize_mapping(mapping): - clean_keys = sanitize_keys([k for k in mapping.keys()]) - for old_key, new_key in clean_keys.items(): - if old_key != new_key: - mapping[new_key] = mapping[old_key] - del mapping[old_key] - - -def clean_all_column_names(adata): - sanitize_df(adata.obs) - sanitize_df(adata.var) - sanitize_mapping(adata.obsm) - - -if __name__ == "__main__": - main() diff --git a/server/converters/h5ad_data_file.py b/server/converters/h5ad_data_file.py new file mode 100644 index 00000000..c8223822 --- /dev/null +++ b/server/converters/h5ad_data_file.py @@ -0,0 +1,250 @@ +import json +import logging +from os import path + +import anndata +import numpy as np +import tiledb + +from server.common.colors import convert_anndata_category_colors_to_cxg_category_colors +from server.common.corpora import corpora_get_props_from_anndata +from server.common.errors import ColorFormatException +from server.common.utils.cxg_constants import CxgConstants +from server.common.utils.cxg_generation_utils import ( + convert_dictionary_to_cxg_group, + convert_dataframe_to_cxg_array, + convert_ndarray_to_cxg_dense_array, + convert_matrix_to_cxg_array, +) +from server.common.utils.matrix_utils import is_matrix_sparse, get_column_shift_encode_for_matrix + + +class H5ADDataFile: + """ Class encapsulating required information about an H5AD datafile that ultimately will be transformed into + another format (currently just CXG is supported). """ + + def __init__( + self, + input_filename, + backed=False, + dataset_title=None, + dataset_about=None, + obs_index_column_name=None, + vars_index_column_name=None, + use_corpora_schema=True, + ): + self.input_filename = input_filename + self.backed = backed + self.dataset_title = dataset_title + self.dataset_about = dataset_about + self.obs_index_column_name = obs_index_column_name + self.vars_index_column_name = vars_index_column_name + + self.use_corpora_schema = use_corpora_schema + + self.validate_input_file_type() + + self.extract_anndata_elements_from_file() + self.extract_metadata_about_dataset() + + self.validate_anndata() + + def to_cxg(self, output_cxg_directory, sparse_threshold, convert_anndata_colors_to_cxg_colors=True): + """ + Writes the following attributes of the anndata to CXG: 1) the metadata as metadata attached to an empty + DenseArray, 2) the obs DataFrame as a DenseArray, 3) the var DataFrame as a DenseArray, 4) all valid + embeddings stored in obsm, each one as a DenseArray, 5) the main X matrix of the anndata as either a + SparseArray or DenseArray based on the `sparse_threshold`, and optionally 6) the column shift of the main X + matrix that might turn an otherwise Dense matrix into a Sparse matrix. + """ + + logging.info("Beginning writing to CXG.") + ctx = tiledb.Ctx( + { + "sm.num_reader_threads": 32, + "sm.num_writer_threads": 32, + "sm.consolidation.buffer_size": 1 * 1024 * 1024 * 1024, + } + ) + + tiledb.group_create(output_cxg_directory, ctx=ctx) + logging.info(f"\t...group created, with name {output_cxg_directory}") + + convert_dictionary_to_cxg_group( + output_cxg_directory, self.generate_cxg_metadata(convert_anndata_colors_to_cxg_colors) + ) + logging.info("\t...dataset metadata saved") + + convert_dataframe_to_cxg_array(output_cxg_directory, "obs", self.obs, self.obs_index_column_name, ctx) + logging.info("\t...dataset obs dataframe saved") + + convert_dataframe_to_cxg_array(output_cxg_directory, "var", self.var, self.var_index_column_name, ctx) + logging.info("\t...dataset var dataframe saved") + + self.write_anndata_embeddings_to_cxg(output_cxg_directory, ctx) + logging.info("\t...dataset embeddings saved") + + self.write_anndata_x_matrix_to_cxg(output_cxg_directory, ctx, sparse_threshold) + logging.info("\t...dataset X matrix saved") + + logging.info("Completed writing to CXG.") + + def write_anndata_x_matrix_to_cxg(self, output_cxg_directory, ctx, sparse_threshold): + matrix_container = f"{output_cxg_directory}/X" + + x_matrix_data = self.anndata.X + is_sparse = is_matrix_sparse(x_matrix_data, sparse_threshold) + if not is_sparse: + col_shift = get_column_shift_encode_for_matrix(x_matrix_data, sparse_threshold) + is_sparse = col_shift is not None + else: + col_shift = None + + if col_shift is not None: + logging.info("Converting matrix X as sparse matrix with column shift encoding") + x_col_shift_name = f"{output_cxg_directory}/X_col_shift" + convert_ndarray_to_cxg_dense_array(x_col_shift_name, col_shift, ctx) + + convert_matrix_to_cxg_array(matrix_container, x_matrix_data, is_sparse, ctx, col_shift) + + tiledb.consolidate(matrix_container, ctx=ctx) + if hasattr(tiledb, "vacuum"): + tiledb.vacuum(matrix_container) + + def write_anndata_embeddings_to_cxg(self, output_cxg_directory, ctx): + def is_valid_embedding(adata, embedding_name, embedding_array): + """ + Returns true if this layout data is a valid array for front-end presentation with the following criteria: + * ndarray, with shape (n_obs, >= 2), dtype float/int/uint + * follows ScanPy embedding naming conventions + * with all values finite or NaN (no +Inf or -Inf) + """ + + is_valid = isinstance(embedding_name, str) and embedding_name.startswith("X_") and len(embedding_name) > 2 + is_valid = is_valid and isinstance(embedding_array, np.ndarray) and embedding_array.dtype.kind in "fiu" + is_valid = is_valid and embedding_array.shape[0] == adata.n_obs and embedding_array.shape[1] >= 2 + is_valid = is_valid and not np.any(np.isinf(embedding_array)) and not np.all(np.isnan(embedding_array)) + return is_valid + + embedding_container = f"{output_cxg_directory}/emb" + tiledb.group_create(embedding_container, ctx=ctx) + + for embedding_name, embedding_values in self.anndata.obsm.items(): + if is_valid_embedding(self.anndata, embedding_name, embedding_values): + embedding_name = f"{embedding_container}/{embedding_name[2:]}" + convert_ndarray_to_cxg_dense_array(embedding_name, embedding_values, ctx) + logging.info(f"\t\t...{embedding_name} embedding created") + + def generate_cxg_metadata(self, convert_anndata_colors_to_cxg_colors): + """ + Return a dictionary containing metadata about CXG dataset. This include data about the version as well as + Corpora schema properties if they exist, among other pieces of metadata. + """ + + cxg_group_metadata = { + "cxg_version": CxgConstants.CXG_VERSION, + "cxg_properties": json.dumps({"title": self.dataset_title, "about": self.dataset_about}), + } + if self.corpora_properties is not None: + cxg_group_metadata["corpora"] = json.dumps(self.corpora_properties) + + if convert_anndata_colors_to_cxg_colors: + try: + cxg_group_metadata["cxg_category_colors"] = json.dumps( + convert_anndata_category_colors_to_cxg_category_colors(self.anndata) + ) + except ColorFormatException: + logging.warning( + "Failed to extract colors from H5AD file! Fix the H5AD file or rerun with " + "--disable-custom-colors. See help for more details." + ) + + return cxg_group_metadata + + def validate_input_file_type(self): + """ + Validate that the input file is of a type that we can handle. Currently the only valid file type is `.h5ad`. + """ + + if not self.input_filename.endswith(".h5ad"): + raise Exception(f"Cannot process input file {self.input_filename}. File must be an H5AD.") + + if self.dataset_title or self.dataset_about: + logging.warning( + "If you convert this dataset into CXG and you explicit specify values for the dataset title metadata " + "or the dataset about metadata, it will override any metadata that is extracted as part of the " + "Corpora schema fields." + ) + + def validate_anndata(self): + if not self.var.index.is_unique: + raise ValueError("Variable index in AnnData object is not unique.") + if not self.obs.index.is_unique: + raise ValueError("Observation index in AnnData object is not unique.") + + def extract_anndata_elements_from_file(self): + logging.info(f"Reading in AnnData dataset: {path.basename(self.input_filename)}") + self.anndata = anndata.read_h5ad(self.input_filename, backed="r" if self.backed else None) + logging.info("Completed reading in AnnData dataset!") + + self.obs = self.transform_dataframe_index_into_column(self.anndata.obs, "obs", self.obs_index_column_name) + self.var = self.transform_dataframe_index_into_column(self.anndata.var, "var", self.vars_index_column_name) + + def extract_metadata_about_dataset(self): + """ + Extract metadata information about the dataset that upon conversion will be saved as group metadata with the + CXG that is generated. This metadata information includes Corpora schema properties, the dataset title and + a link that details more information about the dataset. + """ + + self.corpora_properties = corpora_get_props_from_anndata(self.anndata) if self.use_corpora_schema else None + if self.corpora_properties is None and self.use_corpora_schema: + # If the return value is None, this means that we were not able to figure out what version of the Corpora + # schema the object is using and therefore cannot extract any properties. + raise ValueError("Unknown source file schema version is unsupported.") + + # The title and about properties of the dataset are set by the following order: if they are explicitly defined + # then use the explicit value. If the dataset is a Corpora-schema based schema, then extract the title and about + # from the corpora_properties. Otherwise, use the input filename (only for title, about will be blank). + if self.corpora_properties: + corpora_project_links = self.corpora_properties.get("project_links", []) + corpora_about_link = next( + (link for link in corpora_project_links if (link.get("link_type", None) == "SUMMARY")), {} + ) + else: + corpora_about_link = {} + + filename = path.splitext(path.basename(self.input_filename))[0] + + self.dataset_title = self.dataset_title if self.dataset_title else corpora_about_link.get("link_name", filename) + self.dataset_about = self.dataset_about if self.dataset_about else corpora_about_link.get("link_url") + + def transform_dataframe_index_into_column(self, dataframe, dataframe_name, index_column_name): + """ + Convert the dataframe's index into another column in the dataframe. If an index_column_name is specified, + use that column as the index instead. + """ + + if index_column_name is None: + # Create a unique column name for the index. + suffix = 0 + while f"name_{suffix}" in dataframe.columns: + suffix += 1 + index_column_name = f"name_{suffix}" + + # Turn the index into a normal column + dataframe.rename_axis(index_column_name, inplace=True) + dataframe.reset_index(inplace=True) + + elif index_column_name in dataframe.columns: + # User has specified alternative column for unique names, and it exists + if not dataframe[index_column_name].is_unique: + raise KeyError( + f"Values in {dataframe_name}.{index_column_name} must be unique. Please prepare data to contain " + f"unique values." + ) + else: + raise KeyError(f"Column {index_column_name} does not exist.") + + setattr(self, f"{dataframe_name}_index_column_name", index_column_name) + return dataframe diff --git a/server/test/unit/common/test_api.py b/server/test/unit/common/test_api.py index 3e2c79b9..ee979b1f 100644 --- a/server/test/unit/common/test_api.py +++ b/server/test/unit/common/test_api.py @@ -423,11 +423,7 @@ class EndPointsAnndataAnnotations(unittest.TestCase, EndPointsAnnotations): cls.data, cls.tmp_dir, cls.annotations = data_with_tmp_annotations( MatrixDataType.H5AD, annotations_fixture=True ) - cls._setupClass(cls, [ - "--annotations-file", - cls.annotations.output_file, - cls.data.get_location(), - ]) + cls._setupClass(cls, ["--annotations-file", cls.annotations.output_file, cls.data.get_location(), ]) @classmethod def tearDownClass(cls): diff --git a/server/test/unit/common/test_corpora.py b/server/test/unit/common/test_corpora.py index f502244c..1a0222dc 100644 --- a/server/test/unit/common/test_corpora.py +++ b/server/test/unit/common/test_corpora.py @@ -1,9 +1,10 @@ -import unittest -import anndata import json -import tempfile import shutil +import tempfile +import unittest from http import HTTPStatus + +import anndata import requests from server.common.corpora import ( @@ -104,7 +105,7 @@ class CorporaRESTAPITest(unittest.TestCase): "project_links": json.dumps([ {"link_name": "test link", "link_type": "SUMMARY", "link_url": "https://a.u.r.l/"} ]), - "default_embedding": "X_tsne" + "default_embedding": "X_tsne", } adata.uns.update(corpora_props) adata.write(path) diff --git a/server/test/unit/common/utils/test_cxg_generation_utils.py b/server/test/unit/common/utils/test_cxg_generation_utils.py new file mode 100644 index 00000000..57893913 --- /dev/null +++ b/server/test/unit/common/utils/test_cxg_generation_utils.py @@ -0,0 +1,149 @@ +import json +import unittest +from os import popen, path, mkdir +from shutil import rmtree +from uuid import uuid4 + +import numpy as np +import tiledb +from pandas import Series, DataFrame + +from server.common.utils.cxg_generation_utils import (convert_dictionary_to_cxg_group, convert_dataframe_to_cxg_array, + convert_ndarray_to_cxg_dense_array, convert_matrix_to_cxg_array) + +PROJECT_ROOT = popen("git rev-parse --show-toplevel").read().strip() + + +class TestCxgGenerationUtils(unittest.TestCase): + def setUp(self): + self.testing_cxg_temp_directory = f"{PROJECT_ROOT}/server/test/fixtures/{uuid4()}" + mkdir(self.testing_cxg_temp_directory) + + def tearDown(self): + if path.isdir(self.testing_cxg_temp_directory): + rmtree(self.testing_cxg_temp_directory) + + def test__convert_dictionary_to_cxg_group__writes_successfully(self): + random_dictionary = {"cookies": "chocolate_chip", "brownies": "chocolate", "cake": "double chocolate"} + dictionary_name = "favorite_desserts" + expected_array_directory = f"{self.testing_cxg_temp_directory}/{dictionary_name}" + + convert_dictionary_to_cxg_group(self.testing_cxg_temp_directory, random_dictionary, + group_metadata_name=dictionary_name) + + array = tiledb.open(expected_array_directory) + actual_stored_metadata = dict(array.meta.items()) + + self.assertTrue(path.isdir(expected_array_directory)) + self.assertTrue(isinstance(array, tiledb.DenseArray)) + self.assertEqual(random_dictionary, actual_stored_metadata) + + def test__convert_dataframe_to_cxg_array__writes_successfully(self): + random_int_category = Series(data=[3, 1, 2, 4], dtype=np.int64) + random_bool_category = Series(data=[True, True, False, True], dtype=np.bool_) + random_dataframe_name = f"random_dataframe_{uuid4()}" + random_dataframe = DataFrame(data={"int_category": random_int_category, "bool_category": random_bool_category}) + + convert_dataframe_to_cxg_array(self.testing_cxg_temp_directory, random_dataframe_name, random_dataframe, + "int_category", tiledb.Ctx()) + + expected_array_directory = f"{self.testing_cxg_temp_directory}/{random_dataframe_name}" + expected_array_metadata = { + "cxg_schema": json.dumps({"int_category": {"type": "int32"}, "bool_category": {"type": "boolean"}, + "index": "int_category"})} + + actual_stored_dataframe_array = tiledb.open(expected_array_directory) + actual_stored_dataframe_metadata = dict(actual_stored_dataframe_array.meta.items()) + + self.assertTrue(path.isdir(expected_array_directory)) + self.assertTrue(isinstance(actual_stored_dataframe_array, tiledb.DenseArray)) + self.assertDictEqual(expected_array_metadata, actual_stored_dataframe_metadata) + self.assertTrue((actual_stored_dataframe_array[0:4]["int_category"] == random_int_category.to_numpy()).all()) + self.assertTrue((actual_stored_dataframe_array[0:4]["bool_category"] == random_bool_category.to_numpy()).all()) + + def test__convert_ndarray_to_cxg_dense_array__writes_successfully(self): + ndarray = np.random.rand(3, 2) + ndarray_name = f"{self.testing_cxg_temp_directory}/awesome_ndarray_{uuid4()}" + + convert_ndarray_to_cxg_dense_array(ndarray_name, ndarray, tiledb.Ctx()) + + actual_stored_array = tiledb.open(ndarray_name) + + self.assertTrue(path.isdir(ndarray_name)) + self.assertTrue(isinstance(actual_stored_array, tiledb.DenseArray)) + self.assertTrue((actual_stored_array[:, :] == ndarray).all()) + + def test__convert_matrix_to_cxg_array__dense_array_writes_successfully(self): + matrix = np.float32(np.random.rand(3, 2)) + matrix_name = f"{self.testing_cxg_temp_directory}/awesome_matrix_{uuid4()}" + + convert_matrix_to_cxg_array(matrix_name, matrix, False, tiledb.Ctx()) + + actual_stored_array = tiledb.open(matrix_name) + + self.assertTrue(path.isdir(matrix_name)) + self.assertTrue(isinstance(actual_stored_array, tiledb.DenseArray)) + self.assertTrue((actual_stored_array[:, :] == matrix).all()) + + def test__convert_matrix_to_cxg_array__sparse_array_only_store_nonzeros_empty_array(self): + matrix = np.zeros([3, 2]) + matrix_name = f"{self.testing_cxg_temp_directory}/awesome_zero_matrix_{uuid4()}" + + convert_matrix_to_cxg_array(matrix_name, matrix, True, tiledb.Ctx()) + + actual_stored_array = tiledb.open(matrix_name) + + self.assertTrue(path.isdir(matrix_name)) + self.assertTrue(isinstance(actual_stored_array, tiledb.SparseArray)) + self.assertTrue(actual_stored_array[:, :][''].size == 0) + + def test__convert_matrix_to_cxg_array__sparse_array_only_store_nonzeros(self): + matrix = np.zeros([3, 3]) + matrix[0, 0] = 1 + matrix[1, 1] = 1 + matrix[2, 2] = 2 + matrix_name = f"{self.testing_cxg_temp_directory}/awesome_sparse_matrix_{uuid4()}" + + convert_matrix_to_cxg_array(matrix_name, matrix, True, tiledb.Ctx()) + + actual_stored_array = tiledb.open(matrix_name) + + self.assertTrue(path.isdir(matrix_name)) + self.assertTrue(isinstance(actual_stored_array, tiledb.SparseArray)) + self.assertTrue(actual_stored_array[0, 0][''] == 1) + self.assertTrue(actual_stored_array[1, 1][''] == 1) + self.assertTrue(actual_stored_array[2, 2][''] == 2) + self.assertTrue(actual_stored_array[:, :][''].size == 3) + + def test__convert_matrix_to_cxg_array__sparse_array_with_column_encoding_empty_array(self): + matrix_name = f"{self.testing_cxg_temp_directory}/awesome_column_shift_matrix_{uuid4()}" + matrix = np.ones((3, 2)) + # The column shift will be equal to the matrix since subtracting the column shift from the matrix will create + # a matrix of zeros which is sparse. + column_shift = np.ones((3, 2)) + + convert_matrix_to_cxg_array(matrix_name, matrix, True, tiledb.Ctx(), + column_shift_for_sparse_encoding=column_shift) + + actual_stored_array = tiledb.open(matrix_name) + + self.assertTrue(path.isdir(matrix_name)) + self.assertTrue(isinstance(actual_stored_array, tiledb.SparseArray)) + self.assertTrue(actual_stored_array[:, :][''].size == 0) + + def test__convert_matrix_to_cxg_array__sparse_array_with_column_encoding_partial_array(self): + matrix_name = f"{self.testing_cxg_temp_directory}/awesome_column_shift_matrix_{uuid4()}" + matrix = np.ones((2, 2)) + # Only column shift the first column of ones. + column_shift = np.array([[1, 0], [1, 0]]) + + convert_matrix_to_cxg_array(matrix_name, matrix, True, tiledb.Ctx(), + column_shift_for_sparse_encoding=column_shift) + + actual_stored_array = tiledb.open(matrix_name) + + self.assertTrue(path.isdir(matrix_name)) + self.assertTrue(isinstance(actual_stored_array, tiledb.SparseArray)) + self.assertTrue(actual_stored_array[0, 1][''] == 1) + self.assertTrue(actual_stored_array[1, 1][''] == 1) + self.assertTrue(actual_stored_array[:, :][''].size == 2) diff --git a/server/test/unit/common/utils/test_type_conversion_utils.py b/server/test/unit/common/utils/test_type_conversion_utils.py index a5bb0395..9778bfca 100644 --- a/server/test/unit/common/utils/test_type_conversion_utils.py +++ b/server/test/unit/common/utils/test_type_conversion_utils.py @@ -2,10 +2,10 @@ import unittest from unittest.mock import patch import numpy as np -from pandas import Series +from pandas import Series, DataFrame from server.common.utils.type_conversion_utils import can_cast_to_float32, can_cast_to_int32, get_dtype_of_array, \ - get_schema_type_hint_of_array + get_schema_type_hint_of_array, get_dtypes_and_schemas_of_dataframe class TestTypeConversionUtils(unittest.TestCase): @@ -119,3 +119,17 @@ class TestTypeConversionUtils(unittest.TestCase): i=test_type_index): array = Series(data=[], dtype=types[test_type_index]) self.assertEqual(get_schema_type_hint_of_array(array), expected_schema_hints[test_type_index]) + + def test__get_dtypes_and_schemas_of_dataframe__dtype_and_schema_returns_as_expected(self): + float_array = Series(data=[1, 2, 3], dtype=np.dtype(np.float64)) + category_array = Series(data=["a", "b", "b"], dtype="category") + dataframe = DataFrame({"float_array": float_array, "category_array": category_array}) + + expected_data_types_dict = {"float_array": np.float32, "category_array": np.unicode} + expected_schema_type_hints_dict = {"float_array": {"type": "float32"}, + "category_array": {"type": "categorical", "categories": ["a", "b"]}} + + actual_dataframe_data_types, actual_dataframe_schema_type_hints = get_dtypes_and_schemas_of_dataframe(dataframe) + + self.assertEqual(expected_data_types_dict, actual_dataframe_data_types) + self.assertEqual(expected_schema_type_hints_dict, actual_dataframe_schema_type_hints) diff --git a/server/test/unit/compute/test_diffexp_cxg.py b/server/test/unit/compute/test_diffexp_cxg.py index 7a6aabdd..30dc7472 100644 --- a/server/test/unit/compute/test_diffexp_cxg.py +++ b/server/test/unit/compute/test_diffexp_cxg.py @@ -1,14 +1,16 @@ +import os +import tempfile import unittest -from server.data_common.matrix_loader import MatrixDataLoader -from server.test import PROJECT_ROOT, app_config, FIXTURES_ROOT + +import numpy as np + import server.compute.diffexp_cxg as diffexp_cxg import server.compute.diffexp_generic as diffexp_generic -from server.converters.cxgtool import write_cxg, create_cxg_group_metadata -from server.test.performance.create_test_matrix import create_test_h5ad +from server.converters.h5ad_data_file import H5ADDataFile from server.data_common.fbs.matrix import encode_matrix_fbs, decode_matrix_fbs -import numpy as np -import tempfile -import os +from server.data_common.matrix_loader import MatrixDataLoader +from server.test import PROJECT_ROOT, app_config, FIXTURES_ROOT +from server.test.performance.create_test_matrix import create_test_h5ad class DiffExpTest(unittest.TestCase): @@ -98,21 +100,22 @@ class DiffExpTest(unittest.TestCase): 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_dataset_config=dict(embeddings__names=[])) - adata = adaptor_anndata.data + h5adfile_path = os.path.join(dirname, "sparse.h5ad") + create_test_h5ad(h5adfile_path, 2000, 2000, 10, apply_col_shift) + + h5ad_file_to_convert = H5ADDataFile(h5adfile_path, use_corpora_schema=False) sparsename = os.path.join(dirname, "sparse.cxg") - cxg_group_metadata = create_cxg_group_metadata(adata=adata, basefname="sparse.h5ad", title="sparse",) - write_cxg(adata=adata, container=sparsename, cxg_group_metadata=cxg_group_metadata, sparse_threshold=11) + h5ad_file_to_convert.to_cxg(sparsename, 11, True) + + adaptor_anndata = self.load_dataset(h5adfile_path, extra_dataset_config=dict(embeddings__names=[])) + 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") - cxg_group_metadata = create_cxg_group_metadata(adata=adata, basefname="dense.h5ad", title="dense",) - write_cxg(adata=adata, container=densename, cxg_group_metadata=cxg_group_metadata, sparse_threshold=0) + h5ad_file_to_convert.to_cxg(densename, True, 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") diff --git a/server/test/unit/converters/test_cxgtool.py b/server/test/unit/converters/test_cxgtool.py deleted file mode 100644 index 46e26b48..00000000 --- a/server/test/unit/converters/test_cxgtool.py +++ /dev/null @@ -1,41 +0,0 @@ -import shutil -import unittest - -import anndata - -from server.common.data_locator import DataLocator -from server.converters.cxgtool import write_cxg, create_cxg_group_metadata -from server.data_cxg.cxg_adaptor import CxgAdaptor -from server.test import PROJECT_ROOT, app_config, random_string -from server.test.fixtures.fixtures import pbmc3k_colors - - -class TestCxgAdaptor(unittest.TestCase): - def setUp(self) -> None: - self.fixtures = [] - - def tearDown(self) -> None: - try: - for data_locator in self.fixtures: - print("REMOVING ", data_locator) - shutil.rmtree(data_locator) - except FileNotFoundError: - pass - - def test_cxg_category_colors(self): - data = self.convert_pbmc3k(extract_colors=True) - self.assertEqual(data.get_colors(), pbmc3k_colors) - data = self.convert_pbmc3k(extract_colors=False) - self.assertEqual(data.get_colors(), {}) - - def convert_pbmc3k(self, **kwargs): - rand_str = random_string(8) - data_locator = f"/tmp/test_{rand_str}.cxg" - self.fixtures.append(data_locator) - source_h5ad = anndata.read_h5ad(f"{PROJECT_ROOT}/example-dataset/pbmc3k.h5ad") - cxg_group_metadata = create_cxg_group_metadata( - adata=source_h5ad, basefname="pbmc3k.h5ad", title="pbmc3k", **kwargs - ) - write_cxg(adata=source_h5ad, container=data_locator, cxg_group_metadata=cxg_group_metadata) - config = app_config(data_locator) - return CxgAdaptor(DataLocator(data_locator), config) diff --git a/server/test/unit/converters/test_h5ad_data_file.py b/server/test/unit/converters/test_h5ad_data_file.py new file mode 100644 index 00000000..f8587ebd --- /dev/null +++ b/server/test/unit/converters/test_h5ad_data_file.py @@ -0,0 +1,235 @@ +import json +import unittest +from glob import glob +from os import popen, remove, path +from shutil import rmtree +from uuid import uuid4 + +import anndata +import numpy as np +from pandas import Series, DataFrame + +from server.common.utils.corpora_constants import CorporaConstants +from server.converters.h5ad_data_file import H5ADDataFile + +PROJECT_ROOT = popen("git rev-parse --show-toplevel").read().strip() + + +class TestH5ADDataFile(unittest.TestCase): + + def setUp(self): + self.sample_anndata = self._create_sample_anndata_dataset() + self.sample_h5ad_filename = self._write_anndata_to_file(self.sample_anndata) + + self.sample_output_directory = path.splitext(self.sample_h5ad_filename)[0] + ".cxg" + + def tearDown(self): + if self.sample_h5ad_filename: + remove(self.sample_h5ad_filename) + + if path.isdir(self.sample_output_directory): + rmtree(self.sample_output_directory) + + def test__create_h5ad_data_file__non_h5ad_raises_exception(self): + non_h5ad_filename = "my_fancy_dataset.csv" + + with self.assertRaises(Exception) as exception_context: + H5ADDataFile(non_h5ad_filename) + + self.assertIn("File must be an H5AD", str(exception_context.exception)) + + def test__create_h5ad_data_file__assert_warning_outputted_if_dataset_title_or_about_given(self): + with self.assertLogs(level="WARN") as logger: + H5ADDataFile(self.sample_h5ad_filename, dataset_title="My Awesome Dataset", + dataset_about="http://www.awesomedataset.com", use_corpora_schema=False) + + self.assertIn("will override any metadata that is extracted", logger.output[0]) + + def test__create_h5ad_data_file__reads_anndata_successfully(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False) + + self.assertTrue((h5ad_file.anndata.X == self.sample_anndata.X).all()) + self.assertEqual(h5ad_file.anndata.obs.sort_index(inplace=True), + self.sample_anndata.obs.sort_index(inplace=True)) + self.assertEqual(h5ad_file.anndata.var.sort_index(inplace=True), + self.sample_anndata.var.sort_index(inplace=True)) + + for key in h5ad_file.anndata.obsm.keys(): + self.assertIn(key, self.sample_anndata.obsm.keys()) + self.assertTrue((h5ad_file.anndata.obsm[key] == self.sample_anndata.obsm[key]).all()) + + for key in self.sample_anndata.obsm.keys(): + self.assertIn(key, h5ad_file.anndata.obsm.keys()) + self.assertTrue((h5ad_file.anndata.obsm[key] == self.sample_anndata.obsm[key]).all()) + + def test__create_h5ad_data_file__copies_index_of_obs_and_var_to_column(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False) + + # The automatic name chosen for the index should be "name_0" + self.assertNotIn("name_0", self.sample_anndata.obs.columns) + self.assertIn("name_0", h5ad_file.obs.columns) + + self.assertNotIn("name_0", self.sample_anndata.var.columns) + self.assertIn("name_0", h5ad_file.var.columns) + + def test__create_h5ad_data_file__no_copy_if_obs_and_var_index_names_specified(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False, + obs_index_column_name="float_category", vars_index_column_name="int_category") + + self.assertNotIn("name_0", h5ad_file.obs.columns) + self.assertNotIn("name_0", h5ad_file.var.columns) + + def test__create_h5ad_data_file__obs_and_var_index_names_specified_not_unique_raises_exception(self): + + with self.assertRaises(Exception) as exception_context: + H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False, + obs_index_column_name="float_category", vars_index_column_name="bool_category") + + self.assertIn("Please prepare data to contain unique values", str(exception_context.exception)) + + def test__create_h5ad_data_file__obs_and_var_index_names_specified_doesnt_exist_raises_exception(self): + with self.assertRaises(Exception) as exception_context: + H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False, + obs_index_column_name="unknown_category", vars_index_column_name="i_dont_exist") + + self.assertIn("does not exist", str(exception_context.exception)) + + def test__create_h5ad_data_file__extract_about_and_title_from_dataset(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename) + + self.assertEqual(h5ad_file.dataset_title, "random_link_name") + self.assertEqual(h5ad_file.dataset_about, "www.link.com") + + def test__create_h5ad_data_file__inputted_dataset_title_and_about_overrides_extracted(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename, dataset_about="override_about", + dataset_title="override_title") + + self.assertEqual(h5ad_file.dataset_title, "override_title") + self.assertEqual(h5ad_file.dataset_about, "override_about") + + def test__to_cxg__simple_anndata_no_corpora_and_sparse(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False) + h5ad_file.to_cxg(self.sample_output_directory, 100) + + self._validate_expected_generated_list_of_tiledb_files() + + def test__to_cxg__simple_anndata_with_corpora_and_sparse(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename) + h5ad_file.to_cxg(self.sample_output_directory, 100) + + self._validate_expected_generated_list_of_tiledb_files() + + def test__to_cxg__simple_anndata_no_corpora_and_dense(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename, use_corpora_schema=False) + h5ad_file.to_cxg(self.sample_output_directory, 0) + + self._validate_expected_generated_list_of_tiledb_files() + + def test__to_cxg__simple_anndata_with_corpora_and_dense(self): + h5ad_file = H5ADDataFile(self.sample_h5ad_filename) + h5ad_file.to_cxg(self.sample_output_directory, 0) + + self._validate_expected_generated_list_of_tiledb_files() + + def test__to_cxg__with_sparse_column_encoding(self): + anndata = self._create_sample_anndata_dataset() + anndata.X = np.ones((3, 4)) + sparse_with_column_shift_filename = self._write_anndata_to_file(anndata) + + h5ad_file = H5ADDataFile(sparse_with_column_shift_filename) + h5ad_file.to_cxg(self.sample_output_directory, 50) + + self._validate_expected_generated_list_of_tiledb_files(has_column_encoding=True) + + # Clean up + remove(sparse_with_column_shift_filename) + + def _validate_expected_generated_list_of_tiledb_files(self, has_column_encoding=False): + expected_directories, expected_obs_files, expected_var_files = \ + self._get_expected_generated_list_of_tiledb_files() + + for directory in expected_directories: + self.assertTrue(path.isdir(directory)) + + for obs_file in expected_obs_files: + expected_location_of_obs_file = f"{self.sample_output_directory}/obs/*/{obs_file}" + self.assertTrue(path.isfile(glob(expected_location_of_obs_file)[0])) + + for var_file in expected_var_files: + expected_location_of_var_file = f"{self.sample_output_directory}/var/*/{var_file}" + self.assertTrue(path.isfile(glob(expected_location_of_var_file)[0])) + + if has_column_encoding: + self.assertTrue(path.isdir(f"{self.sample_output_directory}/X_col_shift")) + + def _get_expected_generated_list_of_tiledb_files(self): + + # Expected directories + metadata_directory = f"{self.sample_output_directory}/cxg_group_metadata" + main_x_directory = f"{self.sample_output_directory}/X" + overall_embedding_directory = f"{self.sample_output_directory}/emb" + specific_embedding_directory = f"{self.sample_output_directory}/emb/awesome_embedding" + obs_directory = f"{self.sample_output_directory}/obs" + var_directory = f"{self.sample_output_directory}/var" + + # Obs files + obs_files = [] + obs_files.append("name_0.tdb") + obs_files.append("name_0_var.tdb") + obs_files.append("string_category.tdb") + obs_files.append("string_category_var.tdb") + obs_files.append("float_category.tdb") + + # Var files + var_files = [] + var_files.append("name_0.tdb") + var_files.append("name_0_var.tdb") + var_files.append("bool_category.tdb") + var_files.append("int_category.tdb") + + return [metadata_directory, main_x_directory, overall_embedding_directory, specific_embedding_directory, + obs_directory, var_directory], obs_files, var_files + + def _write_anndata_to_file(self, anndata): + temporary_filename = f"{PROJECT_ROOT}/server/test/fixtures/{uuid4()}.h5ad" + anndata.write(temporary_filename) + + return temporary_filename + + def _create_sample_anndata_dataset(self): + # Create X + X = np.random.rand(3, 4) + + # Create obs + random_string_category = Series(data=["a", "b", "b"], dtype="category") + random_float_category = Series(data=[3.2, 1.1, 2.2], dtype=np.float32) + obs_dataframe = DataFrame( + data={"string_category": random_string_category, "float_category": random_float_category}) + obs = obs_dataframe + + # Create vars + random_int_category = Series(data=[3, 1, 2, 4], dtype=np.int32) + random_bool_category = Series(data=[True, True, False, True], dtype=np.bool_) + var_dataframe = DataFrame(data={"int_category": random_int_category, "bool_category": random_bool_category}) + var = var_dataframe + + # Create embeddings + random_embedding = np.random.rand(3, 2) + obsm = {"X_awesome_embedding": random_embedding} + + # Create uns corpora metadata + uns = {} + for metadata_field in CorporaConstants.REQUIRED_SIMPLE_METADATA_FIELDS: + uns[metadata_field] = "random" + + for metadata_field in CorporaConstants.REQUIRED_JSON_ENCODED_METADATA_FIELD: + uns[metadata_field] = json.dumps({"random_key": "random_value"}) + + # Need to carefully set the corpora schema versions in order for tests to pass. + uns["version"] = {"corpora_schema_version": "1.0.0", "corpora_encoding_version": "0.1.0"} + + # Set project links to be a dictionary + uns["project_links"] = json.dumps( + [{"link_name": "random_link_name", "link_url": "www.link.com", "link_type": "SUMMARY"}]) + + return anndata.AnnData(X=X, obs=obs, var=var, obsm=obsm, uns=uns)