Files
cellxgene/server/app/scanpy_engine/scanpy_engine.py
T
Bruce Martin 3660a6cc27 Experimental - manual annotations (#837)
* icons, partway

* redux for values

* onChange

* cancel

* annotations lifecycle for category names

* copy categorical

* edit category

* add Dataframe.withColsFrom

* render user annotations; default add/delete annotation category

* add label name to actions

* category name edit

* error checking improvements

* change schema field isUserAnnotation to writable

* always have an unassigned label; implement delete label

* implement add new label and edit label name

* label current cell selection

* fix select exact bug in crossfilter

* clean up categorical reducer

* fix tests

* remove debugging printf

* implement subset/reset for user annotations

* undo redo support for user annotations

* remove duplicate button from categories

* add modal

* remove obsolete duplicate annotation reducers

* remove old debugging printf

* connect modal to annotation create and dup

* initial full-stack wiring

* finish up end-to-end wiring

* fix existing unit tests

* fix pytests to match new schema API

* remove debugging printfs

* add label file rotation

* remove obsolete comment

* add fbs encode/decode tests

* add tests for writable annotations

* simplify code

* fix hashing bug with FBS encoding

* lint

* fix smoke tests

* improve error checking in Dataframe.withColsFrom

* add unit test for Dataframe.withColsFrom

* add unit test for Dataframe.columns and Dataframe.renameCol

* fix bug in FBS encode, add better error checks, refactor

* add FBS encode/decode test

* add clarifying comment

* clean up action type names; fix state inconsistency in crossfilter update

* change autosave timer to 2.5sec

* sort categorical metadata render order so it remains consistent

* add temporary autogenerated label for add-new-label operation

* fix hover-over label menu interference with cell highlighting

* remove debugging code

* add missing reducer cases & fix typo

* make dataframe memoize more general purpose

* add dev mode for annos

* fix error on select duplicate

* handle zero occupancy categories

* correctly maintain unclipped AND clipped world

* correctly handle zero length FBS matrix and label files

* ensure all writable categorical schema contains an unassigned category

* handle case where building occupancy stack for category with no members

* dialog for creating label, disable button if duplicate or empty

* visually separate writeable

* edit category

* fix edit category name

* remove debugging code

* fix edit annotation label

* visually define unassigned, change options

* Pull in requirements.txt from `master`

* label currently selected cells

* duplicate label

* lint

* fix pytest merge issues

* rename --label-file to --experimental-label-file

* remove debugging console log

* spelling error fix; fix bug found in PR review.

* lint
2019-09-18 07:33:41 -04:00

568 lines
22 KiB
Python

import warnings
import copy
import threading
import numpy as np
import pandas
from pandas.core.dtypes.dtypes import CategoricalDtype
import anndata
from scipy import sparse
from server.app.driver.driver import CXGDriver
from server.app.util.constants import Axis, DEFAULT_TOP_N, MAX_LAYOUTS
from server.app.util.errors import (
FilterError,
JSONEncodingValueError,
PrepareError,
ScanpyFileError,
DisabledFeatureError,
)
from server.app.util.utils import jsonify_scanpy, requires_data
from server.app.scanpy_engine.diffexp import diffexp_ttest
from server.app.util.fbs.matrix import encode_matrix_fbs, decode_matrix_fbs
from server.app.scanpy_engine.labels import read_labels, write_labels
"""
Sort order for methods
1. Initialize
2. Helper
3. Filter
4. Data & Metadata
5. Computation
"""
class ScanpyEngine(CXGDriver):
def __init__(self, data=None, args={}):
super().__init__(data, args)
# lock used to protect label file write ops
self.label_lock = threading.Lock()
if self.data:
self._validate_and_initialize()
def update(self, data=None, args={}):
super().__init__(data, args)
if self.data:
self._validate_and_initialize()
@staticmethod
def _get_default_config():
return {
"layout": [],
"max_category_items": 100,
"obs_names": None,
"var_names": None,
"diffexp_lfc_cutoff": 0.01,
"label_file": None,
}
@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).
"""
for (ax_name, config_name) in ((Axis.OBS, "obs_names"), (Axis.VAR, "var_names")):
name = self.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.config[config_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)
else:
# user specified a non-existent column name
raise KeyError(
f"Annotation name {name}, specified in --{ax_name}-name does not exist."
)
@staticmethod
def _can_cast_to_float32(ann):
if ann.dtype.kind == "f":
if not np.can_cast(ann.dtype, np.float32):
warnings.warn(
f"Annotation {ann.name} will be converted to 32 bit float and may lose precision."
)
return True
return False
@staticmethod
def _can_cast_to_int32(ann):
if ann.dtype.kind in ["i", "u"]:
if np.can_cast(ann.dtype, np.int32):
return True
ii32 = np.iinfo(np.int32)
if ann.min() >= ii32.min and ann.max() <= ii32.max:
return True
return False
@staticmethod
def _get_col_type(col):
dtype = col.dtype
data_kind = dtype.kind
schema = {}
if ScanpyEngine._can_cast_to_float32(col):
schema["type"] = "float32"
elif ScanpyEngine._can_cast_to_int32(col):
schema["type"] = "int32"
elif dtype == np.bool_:
schema["type"] = "boolean"
elif data_kind == "O" and dtype == "object":
schema["type"] = "string"
elif data_kind == "O" and dtype == "category":
schema["type"] = "categorical"
schema["categories"] = dtype.categories.tolist()
else:
raise TypeError(
f"Annotations of type {dtype} are unsupported by cellxgene."
)
return schema
@requires_data
def _create_schema(self):
self.schema = {
"dataframe": {
"nObs": self.cell_count,
"nVar": self.gene_count,
"type": str(self.data.X.dtype),
},
"annotations": {
"obs": {
"index": self.config["obs_names"],
"columns": []
},
"var": {
"index": self.config["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(self._get_col_type(curr_axis[ann]))
self.schema["annotations"][ax]["columns"].append(ann_schema)
for layout in self.config['layout']:
layout_schema = {
"name": layout,
"type": "float32",
"dims": [f"{layout}_0", f"{layout}_1"]
}
self.schema["layout"]["obs"].append(layout_schema)
@requires_data
def get_schema(self):
schema = self.schema # base schema
# add label obs annotations as needed
if self.labels is not None:
schema = copy.deepcopy(schema)
for col in self.labels.columns:
col_schema = {
"name": col,
"writable": True,
}
col_schema.update(self._get_col_type(self.labels[col]))
schema["annotations"]["obs"]["columns"].append(col_schema)
return 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.
self.data = anndata.read_h5ad(lh)
except ValueError:
raise ScanpyFileError(
"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 ScanpyFileError("Out of memory - file is too large for available memory.")
except Exception as e:
raise ScanpyFileError(
f"{e} - file not found or is inaccessible. File must be an .h5ad object. "
f"Please check your input and try again."
)
if self.config["label_file"]:
try:
self.labels = read_labels(self.config["label_file"])
except Exception as e:
raise ScanpyFileError(
f"Error while loading label file: {e}, File must be in the .csv format, please check "
f"your input and try again."
)
else:
self.labels = None
@requires_data
def _validate_and_initialize(self):
# 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(f"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._default_and_validate_layouts()
self._validate_label_file()
self._create_schema()
@requires_data
def _default_and_validate_layouts(self):
""" function:
a) generate list of default layouts, if not already user specified
b) validate layouts are legal. remove/warn on any that are not
c) cap total list of layouts at global const MAX_LAYOUTS
"""
layouts = self.config['layout']
# handle default
if layouts is None or len(layouts) == 0:
# load default layouts from the data.
layouts = [key[2:] for key in self.data.obsm_keys() if type(key) == str and key.startswith("X_")]
if len(layouts) == 0:
raise PrepareError(f"Unable to find any precomputed layouts within the dataset.")
# 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(f"No valid layout data.")
# cap layouts to MAX_LAYOUTS
self.config['layout'] = valid_layouts[0:MAX_LAYOUTS]
@requires_data
def _is_valid_layout(self, 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
* contains only finite values
"""
is_valid = type(arr) == 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 np.all(np.isfinite(arr))
return is_valid
@requires_data
def _validate_data_types(self):
if sparse.isspmatrix(self.data.X) and not sparse.isspmatrix_csc(self.data.X):
warnings.warn(
f"Scanpy data matrix is sparse, but not a CSC (columnar) matrix. "
f"Performance may be improved by using CSC."
)
if self.data.X.dtype != "float32":
warnings.warn(
f"Scanpy data matrix is in {self.data.X.dtype} format not float32. "
f"Precision may be truncated."
)
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"Scanpy 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.config['max_category_items']:
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"
)
@requires_data
def _validate_label_file(self):
"""
labels is None if disabled, empty if enabled by no data
"""
if self.labels is None or self.labels.empty:
return
# all lables must have a name, which must be unique and not used in obs column names
if not self.labels.columns.is_unique:
raise KeyError(f"All column names specified in {self.config['label_file']} must be unique.")
duplicate_columns = list(set(self.labels.columns) & set(self.data.obs.columns))
if len(duplicate_columns) > 0:
raise KeyError(f"Labels file may not contain column names which overlap "
f"with h5ad obs columns {duplicate_columns}")
# labels must have same count as obs annotations
if self.labels.shape[0] != self.data.obs.shape[0]:
raise ValueError("Labels file must have same number of rows as h5ad file.")
@staticmethod
def _annotation_filter_to_mask(filter, d_axis, count):
mask = np.ones((count,), dtype=bool)
for v in filter:
if d_axis[v["name"]].dtype.name in ["boolean", "category", "object"]:
key_idx = np.in1d(getattr(d_axis, v["name"]), v["values"])
mask = np.logical_and(mask, key_idx)
else:
min_ = v.get("min", None)
max_ = v.get("max", None)
if min_ is not None:
key_idx = (getattr(d_axis, v["name"]) >= min_).ravel()
mask = np.logical_and(mask, key_idx)
if max_ is not None:
key_idx = (getattr(d_axis, v["name"]) <= max_).ravel()
mask = np.logical_and(mask, key_idx)
return mask
@staticmethod
def _index_filter_to_mask(filter, count):
mask = np.zeros((count,), dtype=bool)
for i in filter:
if type(i) == list:
mask[i[0]: i[1]] = True
else:
mask[i] = True
return mask
@staticmethod
def _axis_filter_to_mask(filter, d_axis, count):
mask = np.ones((count,), dtype=bool)
if "index" in filter:
mask = np.logical_and(
mask, ScanpyEngine._index_filter_to_mask(filter["index"], count)
)
if "annotation_value" in filter:
mask = np.logical_and(
mask,
ScanpyEngine._annotation_filter_to_mask(
filter["annotation_value"], d_axis, count
),
)
return mask
@requires_data
def _filter_to_mask(self, filter, use_slices=True):
if use_slices:
obs_selector = slice(0, self.data.n_obs)
var_selector = slice(0, self.data.n_vars)
else:
obs_selector = None
var_selector = None
if filter is not None:
if Axis.OBS in filter:
obs_selector = self._axis_filter_to_mask(
filter["obs"], self.data.obs, self.data.n_obs
)
if Axis.VAR in filter:
var_selector = self._axis_filter_to_mask(
filter["var"], self.data.var, self.data.n_vars
)
return obs_selector, var_selector
@requires_data
def annotation_to_fbs_matrix(self, axis, fields=None):
if axis == Axis.OBS:
if self.labels is not None and not self.labels.empty:
df = pandas.concat([self.data.obs, self.labels], axis=1, join_axes=[self.data.obs.index], copy=False)
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)
@requires_data
def annotation_put_fbs(self, axis, fbs):
fname = self.config["label_file"]
if not fname or self.labels is None:
raise DisabledFeatureError("Writable annotations are not enabled")
if axis != Axis.OBS:
raise ValueError("Only OBS dimension access is supported")
new_label_df = decode_matrix_fbs(fbs)
# if any of the new column labels overlap with our existing labels, raise error
duplicate_columns = list(set(new_label_df.columns) & set(self.data.obs.columns))
if not new_label_df.columns.is_unique or len(duplicate_columns) > 0:
raise KeyError(f"Labels file may not contain column names which overlap "
f"with h5ad obs columns {duplicate_columns}")
# update our internal state and save it. Multi-threading often enabled,
# so treat this as a critical section critical section.
with self.label_lock:
self.labels = new_label_df
write_labels(fname, self.labels)
return jsonify_scanpy({"status": "OK"})
@staticmethod
def slice_columns(X, var_mask):
"""
Slice columns from the matrix X, as specified by the mask
Semantically equivalent to X[:, var_mask], but handles sparse
matrices in a more performant manner.
"""
if var_mask is None: # noop
return X
if sparse.issparse(X): # use tuned getcol/hstack for performance
indices = np.nonzero(var_mask)[0]
cols = [X.getcol(i) for i in indices]
return sparse.hstack(cols, format="csc")
else: # else, just use standard slicing, which is fine for dense arrays
return X[:, var_mask]
@requires_data
def data_frame_to_fbs_matrix(self, filter, axis):
"""
Retrieves data 'X' and returns in a flatbuffer Matrix.
:param filter: filter: dictionary with filter params
:param axis: string obs or var
:return: flatbuffer Matrix
Caveats:
* currently only supports access on VAR axis
* currently only supports filtering on VAR axis
"""
if axis != Axis.VAR:
raise ValueError("Only VAR dimension access is supported")
try:
obs_selector, var_selector = self._filter_to_mask(filter, use_slices=False)
except (KeyError, IndexError, TypeError) as e:
raise FilterError(f"Error parsing filter: {e}") from e
if obs_selector is not None:
raise FilterError("filtering on obs unsupported")
# Currently only handles VAR dimension
X = self.slice_columns(self.data._X, var_selector)
return encode_matrix_fbs(X, col_idx=np.nonzero(var_selector)[0], row_idx=None)
@requires_data
def diffexp_topN(self, obsFilterA, obsFilterB, top_n=None, interactive_limit=None):
if Axis.VAR in obsFilterA or Axis.VAR in obsFilterB:
raise FilterError("Observation filters may not contain vaiable conditions")
try:
obs_mask_A = self._axis_filter_to_mask(
obsFilterA["obs"], self.data.obs, self.data.n_obs
)
obs_mask_B = self._axis_filter_to_mask(
obsFilterB["obs"], self.data.obs, self.data.n_obs
)
except (KeyError, IndexError) as e:
raise FilterError(f"Error parsing filter: {e}") from e
if top_n is None:
top_n = DEFAULT_TOP_N
result = diffexp_ttest(
self.data, obs_mask_A, obs_mask_B, top_n, self.config['diffexp_lfc_cutoff']
)
try:
return jsonify_scanpy(result)
except ValueError:
raise JSONEncodingValueError(
"Error encoding differential expression to JSON"
)
@requires_data
def layout_to_fbs_matrix(self):
"""
Return the default 2-D layout for cells as a FBS Matrix.
Caveats:
* does not support filtering
* only returns Matrix in columnar layout
All embeddings must be individually centered & scaled (isotropically)
to a [0, 1] range.
"""
try:
layout_data = []
for layout in self.config["layout"]:
full_embedding = self.data.obsm[f"X_{layout}"]
embedding = full_embedding[:, :2]
# scale isotropically
min = embedding.min(axis=0)
max = embedding.max(axis=0)
scale = np.amax(max - min)
normalized_layout = (embedding - min) / scale
# translate to center on both axis
translate = 0.5 - ((max - min) / scale / 2)
normalized_layout = normalized_layout + translate
normalized_layout = normalized_layout.astype(dtype=np.float32)
layout_data.append(pandas.DataFrame(normalized_layout, columns=[f"{layout}_0", f"{layout}_1"]))
except ValueError as e:
raise PrepareError(
f"Layout has not been calculated using {self.config['layout']}, "
f"please prepare your datafile and relaunch cellxgene") from e
df = pandas.concat(layout_data, axis=1, copy=False)
return encode_matrix_fbs(df, col_idx=df.columns, row_idx=None)