import warnings import anndata import numpy as np from packaging import version from pandas.core.dtypes.dtypes import CategoricalDtype from scipy import sparse import server.common.compute.diffexp_generic as diffexp_generic import server.common.compute.estimate_distribution as estimate_distribution from server.common.colors import convert_anndata_category_colors_to_cxg_category_colors from server.common.constants import Axis, MAX_LAYOUTS, XApproximateDistribution from server.common.corpora import corpora_get_props_from_anndata from server.common.errors import PrepareError, DatasetAccessError from server.common.utils.type_conversion_utils import get_schema_type_hint_of_array from server.data_common.data_adaptor import DataAdaptor from server.common.fbs.matrix import encode_matrix_fbs anndata_version = version.parse(str(anndata.__version__)).release def anndata_version_is_pre_070(): major = anndata_version[0] minor = anndata_version[1] if len(anndata_version) > 1 else 0 return major == 0 and minor < 7 class AnndataAdaptor(DataAdaptor): def __init__(self, data_locator, app_config=None, dataset_config=None): super().__init__(data_locator, app_config, dataset_config) self.data = None self.X_approximate_distribution = None self._load_data(data_locator) self._validate_and_initialize() def cleanup(self): pass @staticmethod def pre_load_validation(data_locator): if data_locator.islocal(): # if data locator is local, apply file system conventions and other "cheap" # validation checks. If a URI, defer until we actually fetch the data and # try to read it. Many of these tests don't make sense for URIs (eg, extension- # based typing). if not data_locator.exists(): raise DatasetAccessError("does not exist") if not data_locator.isfile(): raise DatasetAccessError("is not a file") @staticmethod def file_size(data_locator): return data_locator.size() if data_locator.islocal() else 0 @staticmethod def open(data_locator, app_config, dataset_config=None): return AnndataAdaptor(data_locator, app_config, dataset_config) def get_corpora_props(self): return corpora_get_props_from_anndata(self.data) def get_name(self): return "cellxgene anndata adaptor version" def get_library_versions(self): return dict(anndata=str(anndata.__version__)) @staticmethod def _create_unique_column_name(df, 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: suffix += 1 return f"{col_name_prefix}{suffix}" def _alias_annotation_names(self): """ The front-end relies on the existance of a unique, human-readable index for obs & var (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 to the front-end via the obs_names and var_names config (which is incorporated into the schema). """ self.original_obs_index = self.data.obs.index for ax_name, var_name in ((Axis.OBS, "obs"), (Axis.VAR, "var")): config_name = f"single_dataset__{var_name}_names" parameter_name = f"{var_name}_names" name = getattr(self.server_config, config_name) df_axis = getattr(self.data, str(ax_name)) if name is None: # Default: create unique names from index if not df_axis.index.is_unique: raise KeyError( f"Values in {ax_name}.index must be unique. " "Please prepare data to contain unique index values, or specify an " "alternative with --{ax_name}-name." ) name = self._create_unique_column_name(df_axis.columns, "name_") self.parameters[parameter_name] = name # reset index to simple range; alias name to point at the # previously specified index. df_axis.rename_axis(name, inplace=True) df_axis.reset_index(inplace=True) elif name in df_axis.columns: # User has specified alternative column for unique names, and it exists if not df_axis[name].is_unique: raise KeyError( f"Values in {ax_name}.{name} must be unique. " "Please prepare data to contain unique values." ) df_axis.reset_index(drop=True, inplace=True) self.parameters[parameter_name] = name else: # user specified a non-existent column name raise KeyError(f"Annotation name {name}, specified in --{ax_name}-name does not exist.") def _create_schema(self): self.schema = { "dataframe": { "nObs": self.cell_count, "nVar": self.gene_count, **get_schema_type_hint_of_array(self.data.X), }, "annotations": { "obs": {"index": self.parameters.get("obs_names"), "columns": []}, "var": {"index": self.parameters.get("var_names"), "columns": []}, }, "layout": {"obs": []}, } for ax in Axis: curr_axis = getattr(self.data, str(ax)) for ann in curr_axis: ann_schema = {"name": ann, "writable": False} ann_schema.update(get_schema_type_hint_of_array(curr_axis[ann])) self.schema["annotations"][ax]["columns"].append(ann_schema) for layout in self.get_embedding_names(): layout_schema = {"name": layout, "type": "float32", "dims": [f"{layout}_0", f"{layout}_1"]} self.schema["layout"]["obs"].append(layout_schema) def get_schema(self): return self.schema def _load_data(self, data_locator): # as of AnnData 0.6.19, backed mode performs initial load fast, but at the # cost of significantly slower access to X data. try: # there is no guarantee data_locator indicates a local file. The AnnData # API will only consume local file objects. If we get a non-local object, # make a copy in tmp, and delete it after we load into memory. with data_locator.local_handle() as lh: # as of AnnData 0.6.19, backed mode performs initial load fast, but at the # cost of significantly slower access to X data. backed = "r" if self.server_config.adaptor__anndata_adaptor__backed else None self.data = anndata.read_h5ad(lh, backed=backed) except ValueError: raise DatasetAccessError( "File must be in the .h5ad format. Please read " "https://github.com/theislab/scanpy_usage/blob/master/170505_seurat/info_h5ad.md to " "learn more about this format. You may be able to convert your file into this format " "using `cellxgene prepare`, please run `cellxgene prepare --help` for more " "information." ) except MemoryError: raise DatasetAccessError("Out of memory - file is too large for available memory.") except Exception: import traceback message = ( "File not found or is inaccessible. File must be an .h5ad object. " "Please check your input and try again." ) if self.server_config.app__verbose: message += f"\n{traceback.format_exc()}" raise DatasetAccessError(message) def _validate_and_initialize(self): if anndata_version_is_pre_070(): warnings.warn( "Use of anndata versions older than 0.7 will have serious issues. Please update to at " "least anndata 0.7 or later." ) # var and obs column names must be unique if not self.data.obs.columns.is_unique or not self.data.var.columns.is_unique: raise KeyError("All annotation column names must be unique.") self._alias_annotation_names() self._validate_data_types() self.cell_count = self.data.shape[0] self.gene_count = self.data.shape[1] self._create_schema() if self.dataset_config.X_approximate_distribution == "auto": """Lazy evaluate the heuristic if we are backed.""" if not self.data.isbacked: self.X_approximate_distribution = estimate_distribution.estimate_approximate_distribution(self.data.X) else: self.X_approximate_distribution = self.dataset_config.X_approximate_distribution # heuristic n_values = self.data.shape[0] * self.data.shape[1] if (n_values > 1e8 and self.server_config.adaptor__anndata_adaptor__backed is True) or (n_values > 5e8): self.parameters.update({"diffexp-may-be-slow": True}) def _is_valid_layout(self, arr): """return True if this layout data is a valid array for front-end presentation: * ndarray, dtype float/int/uint * with shape (n_obs, >= 2) * with all values finite or NaN (no +Inf or -Inf) """ is_valid = type(arr) is np.ndarray and arr.dtype.kind in "fiu" is_valid = is_valid and arr.shape[0] == self.data.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 _validate_data_types(self): # The backed API does not support interrogation of the underlying sparsity or sparse matrix type # Fake it by asking for a small subarray and testing it. NOTE: if the user has ignored our # anndata <= 0.7 warning, opted for the --backed option, and specified a large, sparse dataset, # this "small" indexing request will load the entire X array. This is due to a bug in anndata<=0.7 # which will load the entire X matrix to fullfill any slicing request if X is sparse. See # user warning in _load_data(). X0 = self.data.X[0, 0:1] if sparse.isspmatrix(X0) and not sparse.isspmatrix_csc(X0): warnings.warn( "Anndata data matrix is sparse, but not a CSC (columnar) matrix. " "Performance may be improved by using CSC." ) if self.data.X.dtype > np.dtype(np.float32): warnings.warn( f"Anndata data matrix is in {self.data.X.dtype} format not float32. " f"Precision may be truncated." ) if self.data.X.dtype < np.float32: if self.data.isbacked: raise DatasetAccessError( f"Data matrix in {self.data.X.dtype} format is not supported in backed mode." " Please reload without --backed, or convert matrix to float32" ) warnings.warn( f"Anndata data matrix is in unsupported {self.data.X.dtype} format -- will be cast to float32" ) self.data.X = self.data.X.astype(np.float32) for ax in Axis: curr_axis = getattr(self.data, str(ax)) for ann in curr_axis: datatype = curr_axis[ann].dtype downcast_map = { "int64": "int32", "uint32": "int32", "uint64": "int32", "float64": "float32", } if datatype in downcast_map: warnings.warn( f"Anndata annotation {ax}:{ann} is in unsupported format: {datatype}. " f"Data will be downcast to {downcast_map[datatype]}." ) if isinstance(datatype, CategoricalDtype): category_num = len(curr_axis[ann].dtype.categories) if category_num > 500 and category_num > self.dataset_config.presentation__max_categories: warnings.warn( f"{str(ax).title()} annotation '{ann}' has {category_num} categories, this may be " f"cumbersome or slow to display. We recommend setting the " f"--max-category-items option to 500, this will hide categorical " f"annotations with more than 500 categories in the UI" ) def annotation_to_fbs_matrix(self, axis, fields=None, labels=None): if axis == Axis.OBS: if labels is not None and not labels.empty: df = self.data.obs.join(labels, self.parameters.get("obs_names")) else: df = self.data.obs else: df = self.data.var if fields is not None and len(fields) > 0: df = df[fields] return encode_matrix_fbs(df, col_idx=df.columns) def get_embedding_names(self): """ Return pre-computed embeddings. function: a) generate list of default layouts b) validate layouts are legal. remove/warn on any that are not c) cap total list of layouts at global const MAX_LAYOUTS """ # load default layouts from the data. layouts = self.dataset_config.embeddings__names if layouts is None or len(layouts) == 0: layouts = [key[2:] for key in self.data.obsm_keys() if type(key) is str and key.startswith("X_")] # remove invalid layouts valid_layouts = [] obsm_keys = self.data.obsm_keys() for layout in layouts: layout_name = f"X_{layout}" if layout_name not in obsm_keys: warnings.warn(f"Ignoring unknown layout name: {layout}.") elif not self._is_valid_layout(self.data.obsm[layout_name]): warnings.warn(f"Ignoring layout due to malformed shape or data type: {layout}") else: valid_layouts.append(layout) if len(valid_layouts) == 0: raise PrepareError("No valid layout data.") # cap layouts to MAX_LAYOUTS return valid_layouts[0:MAX_LAYOUTS] def get_embedding_array(self, ename, dims=2): full_embedding = self.data.obsm[f"X_{ename}"] return full_embedding[:, 0:dims] def compute_diffexp_ttest(self, maskA, maskB, top_n=None, lfc_cutoff=None): if top_n is None: top_n = self.dataset_config.diffexp__top_n if lfc_cutoff is None: lfc_cutoff = self.dataset_config.diffexp__lfc_cutoff return diffexp_generic.diffexp_ttest(self, maskA, maskB, top_n, lfc_cutoff) def get_colors(self): return convert_anndata_category_colors_to_cxg_category_colors(self.data) def get_X_array(self, obs_mask=None, var_mask=None): # H5Py does not support boolean indexing (masks), so convert to integer indexing # when backed (ie, when AnnData is using H5Py indexing) if obs_mask is None: obs_mask = slice(None) elif self.data.isbacked and obs_mask.dtype == bool: obs_mask = obs_mask.nonzero()[0] if var_mask is None: var_mask = slice(None) elif self.data.isbacked and var_mask.dtype == bool: var_mask = var_mask.nonzero()[0] X = self.data.X[obs_mask, var_mask] return X def get_X_approximate_distribution(self) -> XApproximateDistribution: """return the approximate distribution of the X matrix.""" if self.X_approximate_distribution is None: """Not yet evaluated.""" assert self.dataset_config.X_approximate_distribution == "auto" self.data = self.data.to_memory() # loads data self.X_approximate_distribution = estimate_distribution.estimate_approximate_distribution(self.data.X) return self.X_approximate_distribution def get_shape(self): return self.data.shape def query_var_array(self, term_name): return getattr(self.data.var, term_name) def query_obs_array(self, term_name): return getattr(self.data.obs, term_name) def get_obs_index(self): name = self.server_config.single_dataset__obs_names if name is None: return self.original_obs_index else: return self.data.obs[name] def get_obs_columns(self): return self.data.obs.columns def get_obs_keys(self): # return list of keys return self.data.obs.keys().to_list() def get_var_keys(self): # return list of keys return self.data.var.keys().to_list()