diff --git a/client/src/actions/index.js b/client/src/actions/index.js index e3da5b62..f8b2d554 100644 --- a/client/src/actions/index.js +++ b/client/src/actions/index.js @@ -90,12 +90,16 @@ async function _doRequestExpressionData(dispatch, getState, genes) { const state = getState(); const { universe } = state.controls; /* preload data already in cache */ - let expressionData = _.transform(genes, (expData, g) => { - const data = kvCache.get(universe.varDataCache, g); - if (data) { - expData[g] = data; - } - }); // --> { gene: data } + let expressionData = _.transform( + genes, + (expData, g) => { + const data = kvCache.get(universe.varDataCache, g); + if (data) { + expData[g] = data; + } + }, + {} + ); // --> { gene: data } /* make a list of genes for which we do not have data */ const genesToFetch = _.filter(genes, g => expressionData[g] === undefined); @@ -119,7 +123,6 @@ async function _doRequestExpressionData(dispatch, getState, genes) { }), headers: new Headers({ accept: "application/json", - "Accept-Encoding": "gzip, deflate, br", "Content-Type": "application/json" }) } @@ -239,7 +242,6 @@ const requestDifferentialExpression = (set1, set2, num_genes = 10) => async ( method: "POST", headers: new Headers({ Accept: "application/json", - "Accept-Encoding": "gzip, deflate, br", "Content-Type": "application/json" }), body: JSON.stringify({ diff --git a/client/src/components/brushableHistogram/index.js b/client/src/components/brushableHistogram/index.js index 838c38d4..b212a16e 100644 --- a/client/src/components/brushableHistogram/index.js +++ b/client/src/components/brushableHistogram/index.js @@ -248,9 +248,9 @@ class HistogramBrush extends React.Component { colorAccessor, isUserDefined, isDiffExp, - avgDiff, - set1AvgExp, - set2AvgExp, + logFoldChange, + pval, + pvalAdj, scatterplotXXaccessor, scatterplotYYaccessor, zebra @@ -337,25 +337,17 @@ class HistogramBrush extends React.Component { }} > - 1: - {` ${set1AvgExp.toPrecision(2)}`} + log fold change: + {` ${logFoldChange.toPrecision(4)}`} - 2: - {` ${set2AvgExp.toPrecision(2)}`} - - - {`Av. Diff: ${avgDiff.toFixed(2)}`} + p-value (adj): + {pvalAdj < 0.0001 ? " < 0.0001" : ` ${pvalAdj.toFixed(4)}`} ) : null} diff --git a/client/src/components/geneExpression/index.js b/client/src/components/geneExpression/index.js index 04f6b871..856d3ddb 100644 --- a/client/src/components/geneExpression/index.js +++ b/client/src/components/geneExpression/index.js @@ -152,9 +152,9 @@ class GeneExpression extends React.Component { zebra={index % 2 === 0} ranges={d3.extent(values)} isDiffExp - avgDiff={value[1]} - set1AvgExp={value[4]} - set2AvgExp={value[5]} + logFoldChange={value[1]} + pval={value[2]} + pvalAdj={value[3]} /> ); }) diff --git a/client/src/util/actionHelpers.js b/client/src/util/actionHelpers.js index 7ee68c59..7ef3c535 100644 --- a/client/src/util/actionHelpers.js +++ b/client/src/util/actionHelpers.js @@ -31,8 +31,7 @@ export const doJsonRequest = async url => { const res = await fetch(url, { method: "get", headers: new Headers({ - "Content-Type": "application/json", - "Accept-Encoding": "gzip, deflate, br" + "Content-Type": "application/json" }) }); if (res.ok && res.headers.get("Content-Type") === "application/json") { diff --git a/docs/REST_API.md b/docs/REST_API.md index 3fb14de6..9ebc71b6 100644 --- a/docs/REST_API.md +++ b/docs/REST_API.md @@ -75,7 +75,7 @@ For a GET URL query parameter: - Annotation name is encoded as `obs:name` or `var:name`[2](#endnote-2). - Enumerated values (string, categorical, boolean) are encoded as option lists, ie, `var:tissue=lung, obs:tumor=true` - Scalar values (int32, float32) are encoded as ranges, ie, `obs:num_reads=1000,10000` where either min or max may be replaced with an asterisk to indicate a half-open range. -- Index filters are not be allowed within GET URL query parameter filters +- Index filters are not allowed within GET URL query parameter filters - Logically, filters are ANDed, except for repeated annotation names which are ORed. For example, `?X=A&X=B&Y=1` is evaluated as `((X==A or X==B) and Y==1)` Example selection for _lung_ and _heart_ tissue with more than 1000 reads: @@ -505,10 +505,10 @@ Generate differential expression (DE) statistics for two specified subsets of da Two modes are provided: -- Return top N differentially expressed variables (genes) -- Return DE for caller-provided variable filter (future) +- `topN`: return top N differentially expressed variables (across all variables) +- `varFilter`: return DE for caller-provided variable filter (_future_) -Both modes perform calculations using a subset of observations, where each subset is defined by an observation filter (`set1` and `set2`). +Both modes perform calculations using a subset of observations, where each subset is defined by an observation filter (`set1` and `set2`). These filters must not include a variable filter. If differential expression is not supported by the server, must return an HTTP 501 response. If, in the view of the server, the request will exceed a reasonable interactive time period, must immediately return HTTP 403 error (error return _before_ attempting computation). @@ -568,24 +568,22 @@ If differential expression is not supported by the server, must return an HTTP 5 **Response body:** -- For 200 Success, differential expression statistics returned as array of arrays sorted by obs index, where each contains the following values: +- For 200 Success, differential expression statistics returned as array of arrays sorted by varindex, where each contains the following values: - **varIndex**: variable index for the computed results - - **avgDiff**: log fold-change of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group, + - **logfoldchange**: log fold-change of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group, - **pVal**: unadjusted p-value, - - **pValAdj**: Adjusted p-value, based on bonferroni correction using all genes in the original dataset), - - **set1AvgExp:** average expression value for all observations in set 1, - - **set2AvgExp**: average expression value for all observations in set 2 + - **pValAdj**: adjusted p-value Statistics are encoded as an array of arrays, with fields ordered as: - _varIndex_, _avgDiff_, _pVal_, _pValAdj_, _set1AvgExp_, _set2AvgExp_ + _varIndex_, _logfoldchange_, _pVal_, _pValAdj_ For example: ``` [ - [ 328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9 ], + [ 1720, 2.4679039, 2.3124478092035228e-175, 4.250279073316075e-172 ] // ... ] ``` @@ -616,8 +614,8 @@ POST /diffexp/obs 200 - Success { "diffexp": [ - [ 328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9 ], - // [ varIdx, avgDiff, pVal, pValAdj, set1AvgExp, set2AvgExp ], + [ 328, -2.569489, 2.655706e-63, 3.642036e-57 ], + // [ varIdx, logfoldchange, pVal, pValAdj ], // ... ] } @@ -690,10 +688,9 @@ Routes: - `GET /schema` - `GET /annotations/obs` - `GET /annotations/var` -- `GET /layout/obs` +- `GET /layout/obs` - get the default layout - `PUT /data/obs` - request will contain a filter by var `name` -- `POST /diffexp/obs` - mode `topN`, typically with a couple of 10, and two sets defined by an obs index filter (`{ filter: { obs: { index: [...] } } }`) -- `PUT /layout/obs` - (_coming soon_) request will contain a filter by obs index +- `POST /diffexp/obs` - mode `topN`, typically with a `count` of 10, and two sets defined by an obs index filter (`{ filter: { obs: { index: [...] } } }`) Requests include the following content negotiation headers: diff --git a/server/app/driver/driver.py b/server/app/driver/driver.py index 3fe0e6ca..b4693336 100644 --- a/server/app/driver/driver.py +++ b/server/app/driver/driver.py @@ -83,16 +83,18 @@ class CXGDriver(metaclass=ABCMeta): pass @abstractmethod - def diffexp(self, filter1, filter2, top_n=None, interactive_limit=None): + def diffexp_topN(self, obsFilter1, obsFilter2, top_n=None, interactive_limit=None): """ - Computes the top differentially expressed variables between two observation sets. If dataframes + Computes the top N differentially expressed variables between two observation sets. If mode + is "TOP_N", then stats for the top N + dataframes contain a subset of variables, then statistics for all variables will be returned, otherwise only the top N vars will be returned. - :param filter1: filter: dictionary with filter params for first set of observations - :param filter2: filter: dictionary with filter params for second set of observations + :param obsFilter1: filter: dictionary with filter params for first set of observations + :param obsFilter2: filter: dictionary with filter params for second set of observations :param top_n: Limit results to top N (Top var mode only) :param interactive_limit: -- don't compute if total # genes in dataframes are larger than this - :return: top genes, stats and expression values for variables + :return: top N genes and corresponding stats """ pass diff --git a/server/app/rest_api/rest.py b/server/app/rest_api/rest.py index 7a21b0de..fecd7dd6 100644 --- a/server/app/rest_api/rest.py +++ b/server/app/rest_api/rest.py @@ -555,11 +555,11 @@ class DiffExpObsAPI(Resource): "responses": { "200": { "description": "Statistics are encoded as an array of arrays, with fields ordered as: " - "varIndex, avgDiff, pVal, pValAdj, set1AvgExp, set2AvgExp", + "varIndex, logfoldchange, pVal, pValAdj", "examples": { "application/json": [ - [328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9], - [1250, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9], + [328, -2.569489, 2.655706e-63, 3.642036e-57], + [1250, -2.569489, 2.655706e-63, 3.642036e-57], ] } }, @@ -584,11 +584,12 @@ class DiffExpObsAPI(Resource): except ValueError: return make_response(f"Error: invalid mode option {args['mode']}", HTTPStatus.BAD_REQUEST) # Validate filters - if mode == DiffExpMode.VAR_FILTER: - if "varFilter" not in args: - return make_response("varFilter is required when mode is set to varFilter ", HTTPStatus.BAD_REQUEST) - if Axis.OBS in args["varFilter"]["filter"]: - return make_response("Obs filter not allowed in varFilter", HTTPStatus.BAD_REQUEST) + if mode == DiffExpMode.VAR_FILTER or "varFilter" in args: + # not NOT_IMPLEMENTED + return make_response("mode=varfilter not implemented", HTTPStatus.NOT_IMPLEMENTED) + if mode == DiffExpMode.TOP_N and "count" not in args: + return make_response("mode=topN requires a count parameter", HTTPStatus.BAD_REQUEST) + if "set1" not in args: return make_response("set1 is required.", HTTPStatus.BAD_REQUEST) if Axis.VAR in args["set1"]["filter"]: @@ -598,16 +599,17 @@ class DiffExpObsAPI(Resource): return make_response("Set2 as inverse of set1 is not implemented", HTTPStatus.NOT_IMPLEMENTED) if Axis.VAR in args["set2"]["filter"]: return make_response("Var filter not allowed for set2", HTTPStatus.BAD_REQUEST) + set1_filter = args["set1"]["filter"] set2_filter = args.get("set2", {"filter": {}})["filter"] - if "varFilter" in args: - set1_filter[Axis.VAR] = args["varFilter"]["filter"][Axis.VAR] - set2_filter[Axis.VAR] = args["varFilter"]["filter"][Axis.VAR] - # mode + + # TODO: implement varfilter mode + + # mode=topN count = args.get("count", None) try: - diffexp = current_app.data.diffexp(set1_filter, set2_filter, count, - current_app.data.features["diffexp"]["interactiveLimit"]) + diffexp = current_app.data.diffexp_topN(set1_filter, set2_filter, count, + current_app.data.features["diffexp"]["interactiveLimit"]) except (ValueError, FilterError) as e: return make_response(e.message, HTTPStatus.BAD_REQUEST) except InteractiveError: diff --git a/server/app/scanpy_engine/diffexp.py b/server/app/scanpy_engine/diffexp.py new file mode 100644 index 00000000..d9c0e2a4 --- /dev/null +++ b/server/app/scanpy_engine/diffexp.py @@ -0,0 +1,83 @@ + +import numpy as np +from scipy import sparse, stats + + +# Convenience function which handles sparse data +def _mean_var_n(X): + """ + Two-pass variance calculation. Numerically (more) stable + than naive methods (and same method used by numpy.var()) + https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass + """ + n = X.shape[0] + if sparse.issparse(X): + mean = X.mean(axis=0).A1 + dfm = X - mean + sumsq = np.sum(np.multiply(dfm, dfm), axis=0).A1 + v = sumsq / (n - 1) + else: + mean = X.mean(axis=0) + dfm = X - mean + sumsq = np.sum(np.multiply(dfm, dfm), axis=0) + v = sumsq / (n - 1) + + return mean, v, n + + +def diffexp_ttest(adata, maskA, maskB, top_n=8): + """ + Return differential expression statistics for top N variables, sorted by + t statistic. Implemented as a unequal variance t-test. + + :param adata: anndata dataframe + :param maskA: observation selection mask for set 1 + :param maskB: observation selection mask for set 2 + :param top_n: number of variables to return stats for + :return: for top N genes, [ varindex, logfoldchange, pval, pval_adj ] + """ + + # mean, variance, N + meanA, vA, nA = _mean_var_n(adata._X[maskA]) + meanB, vB, nB = _mean_var_n(adata._X[maskB]) + + # variance / N + vnA = vA / nA + vnB = vB / nB + sum_vn = vnA + vnB + + # degrees of freedom for Welch's t-test + with np.errstate(divide='ignore', invalid='ignore'): + dof = sum_vn**2 / (vnA**2 / (nA - 1) + vnB**2 / (nB - 1)) + dof[np.isnan(dof)] = 1 + + # Welch's t-test score calculation + with np.errstate(divide='ignore', invalid='ignore'): + tscores = (meanA - meanB) / np.sqrt(sum_vn) + tscores[np.isnan(tscores)] = 0 + + # p-value + pvals = stats.t.sf(np.abs(tscores), dof) * 2 + pvals_adj = pvals * adata._X.shape[1] + + # logfoldchanges: log2(meanA / meanB) + logfoldchanges = np.log2(np.abs((meanA + 1e-9) / (meanB + 1e-9))) + + # top n sort + stats_to_sort = np.abs(tscores) + partition = np.argpartition(stats_to_sort, -top_n)[-top_n:] + rel_sort_order = np.argsort(stats_to_sort[partition])[::-1] + vars_indices = np.arange(adata.n_vars, dtype=int) + sort_order = vars_indices[partition][rel_sort_order] + + # top n slice + logfoldchanges_top_n = logfoldchanges[sort_order] + pvals_top_n = pvals[sort_order] + pvals_adj_top_n = pvals_adj[sort_order] + + # varIndex, logfoldchange, pval, pval_adj + result = [[sort_order[i], + logfoldchanges_top_n[i], + pvals_top_n[i], + pvals_adj_top_n[i]] for i in range(top_n)] + return result diff --git a/server/app/scanpy_engine/scanpy_engine.py b/server/app/scanpy_engine/scanpy_engine.py index 788181b2..79dbe328 100644 --- a/server/app/scanpy_engine/scanpy_engine.py +++ b/server/app/scanpy_engine/scanpy_engine.py @@ -4,11 +4,12 @@ import numpy as np from pandas import DataFrame from pandas.core.dtypes.dtypes import CategoricalDtype import scanpy.api as sc -from scipy import stats, sparse +from scipy import sparse from server.app.driver.driver import CXGDriver -from server.app.util.constants import Axis, DEFAULT_TOP_N, DiffExpMode +from server.app.util.constants import Axis, DEFAULT_TOP_N from server.app.util.errors import FilterError, InteractiveError, PrepareError, ScanpyFileError +from server.app.scanpy_engine.diffexp import diffexp_ttest """ Sort order for methods @@ -114,35 +115,6 @@ class ScanpyEngine(CXGDriver): f"that your input and try again.") return result - @staticmethod - def _top_sort(values, sort_order, top_n=None): - """ - Sorts an iterable in sort order limited by top_n - :param values: iterable of values to sort - :param sort_order: ndarray order to sort in - :param top_n: cutoff number to return - :return: values sorted by sort_order limited by top_n - """ - return values[sort_order][:top_n] - - @staticmethod - def _nan_to_one(values): - """ - Replaces NaN values with 1 - :param values: numpy ndarray - :return: ndarray - """ - return np.where(np.isnan(values), 1, values) - - @staticmethod - def _nan_to_zero(values): - """ - Replaces NaN values with 0 - :param values: numpy ndarray - :return: ndarray - """ - return np.where(np.isnan(values), 0, values) - def _validate_data_types(self): if self.data.X.dtype != "float32": warnings.warn(f"Scanpy data matrix is in {self.data.X.dtype} format not float32. " @@ -180,7 +152,7 @@ class ScanpyEngine(CXGDriver): f"`cellxgene prepare --layout {self.layout_method} ` " f"to solve this problem. ") - def filter_dataframe(self, filter, include_uns=False): + def filter_dataframe(self, filter): """ Filter cells from data and return a subset of the data. They can operate on both obs and var dimension with indexing and filtering by annotation value. Filters are combined with the and operator. @@ -189,70 +161,68 @@ class ScanpyEngine(CXGDriver): https://docs.google.com/document/d/1Fxjp1SKtCk7l8QP9-7KAjGXL0eldi_qEnNT0NmlGzXI/edit#heading=h.8qc9q57amldx :param filter: dictionary with filter params - :param include_uns: bool, include unstructured annotations :return: View into scanpy object with cells/genes filtered """ if not filter: return self.data - cells_idx = np.ones((self.cell_count,), dtype=bool) - genes_idx = np.ones((self.gene_count,), dtype=bool) - if Axis.OBS in filter: - if "index" in filter["obs"]: - cells_idx = self._filter_index(filter["obs"]["index"], cells_idx, Axis.OBS) - if "annotation_value" in filter["obs"]: - cells_idx = self._filter_annotation(filter["obs"]["annotation_value"], cells_idx, Axis.OBS) - if Axis.VAR in filter: - if "index" in filter["var"]: - genes_idx = self._filter_index(filter["var"]["index"], genes_idx, Axis.VAR) - if "annotation_value" in filter["var"]: - genes_idx = self._filter_annotation(filter["var"]["annotation_value"], genes_idx, Axis.VAR) - - data = self._slice(self.data, cells_idx, genes_idx) + obs_selector, var_selector = self._filter_to_mask(filter, use_slices=False) + data = self._slice(self.data, obs_selector, var_selector) return data - def _filter_index(self, filter, index, axis): - """ - Filter data based on index. ex. [1, 3, [111:200]] - :param filter: subset of filter dict for obs/var:index - :param index: np logical vector containing true for passing false for failing filter - :param axis: string obs or var - :return: np logical vector for whether the data passes the filter - """ - if axis == Axis.OBS: - count_ = self.cell_count - elif axis == Axis.VAR: - count_ = self.gene_count - idx_filter = np.zeros((count_,), dtype=bool) - for i in filter: - if type(i) == list: - idx_filter[i[0]:i[1]] = True - else: - idx_filter[i] = True - return np.logical_and(index, idx_filter) - - def _filter_annotation(self, filter, index, axis): - """ - Filter data based on annotation value - :param filter: subset of filter dict for obs/var:annotation_value - :param index: np logical vector containing true for passing false for failing filter - :param axis: string obs or var - :return: np logical vector for whether the data passes the filter - """ - d_axis = getattr(self.data, axis.value) + @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"]) - index = np.logical_and(index, key_idx) + 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() - index = np.logical_and(index, key_idx) + mask = np.logical_and(mask, key_idx) if max_ is not None: key_idx = (getattr(d_axis, v["name"]) <= max_).ravel() - index = np.logical_and(index, key_idx) - return index + 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 + + 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 @staticmethod def _slice(data, obs_selector=None, vars_selector=None): @@ -293,16 +263,25 @@ class ScanpyEngine(CXGDriver): [observation ids, val1, val2...] """ try: - df = self.filter_dataframe(filter) - except KeyError as e: + obs_selector, var_selector = self._filter_to_mask(filter) + except (KeyError, IndexError) as e: raise FilterError(f"Error parsing filter: {e}") from e - df_axis = getattr(df, axis) - if not fields: - fields = df_axis.columns.tolist() - result = { - "names": fields, - "data": DataFrame(df_axis[fields]).to_records(index=True).tolist() - } + if axis == Axis.OBS: + obs = self.data.obs[obs_selector] + if not fields: + fields = obs.columns.tolist() + result = { + "names": fields, + "data": DataFrame(obs[fields]).to_records(index=True).tolist() + } + else: + var = self.data.var[var_selector] + if not fields: + fields = var.columns.tolist() + result = { + "names": fields, + "data": DataFrame(var[fields]).to_records(index=True).tolist() + } return result def data_frame(self, filter, axis): @@ -316,85 +295,38 @@ class ScanpyEngine(CXGDriver): } """ try: - slice = self.filter_dataframe(filter) - except KeyError as e: + obs_selector, var_selector = self._filter_to_mask(filter) + except (KeyError, IndexError) as e: raise FilterError(f"Error parsing filter: {e}") from e - # convert sparse slice to dense - X = slice._X.toarray() if sparse.issparse(slice._X) else slice._X + _X = self.data._X[obs_selector, var_selector] + if sparse.issparse(_X): + _X = _X.toarray() + var_index_sliced = self.data.var.index[var_selector] + obs_index_sliced = self.data.obs.index[obs_selector] if axis == Axis.OBS: result = { - "var": slice.var.index.tolist(), - "obs": DataFrame(X, index=slice.obs.index).to_records(index=True).tolist() + "var": var_index_sliced.tolist(), + "obs": DataFrame(_X, index=obs_index_sliced).to_records(index=True).tolist() } else: result = { - "obs": slice.obs.index.tolist(), - "var": DataFrame(X.T, index=slice.var.index).to_records(index=True).tolist() + "obs": obs_index_sliced.tolist(), + "var": DataFrame(_X.T, index=var_index_sliced).to_records(index=True).tolist() } return result - def diffexp(self, filter1, filter2, top_n=None, interactive_limit=None): - """ - Computes the top differentially expressed variables between two observation sets. If dataframes - contain a subset of variables, then statistics for all variables will be returned, otherwise - only the top N vars will be returned. - :param filter1: filter: dictionary with filter params for first set of observations - :param filter2: filter: dictionary with filter params for second set of observations - :param top_n: Limit results to top N (Top var mode only) - :param interactive_limit: -- don't compute if total # genes in dataframes are larger than this - :return: top genes, stats and expression values for variables - """ + 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: - df1 = self.filter_dataframe(filter1) - except KeyError as e: - raise FilterError(f"Error parsing filter for set 1: {e}") from e - # TODO df2 should be inverse if not filter2 provided - try: - df2 = self.filter_dataframe(filter2) - except KeyError as e: - raise FilterError(f"Error parsing filter for set 2: {e}") from e - # If not the same genes, test is wrong! - if np.any(df1.var.index != df2.var.index): - raise ValueError("Variables ares not the same in set1 and set2") - if interactive_limit and df1.shape[0] + df2.shape[0] > interactive_limit: - raise InteractiveError("Size of set 1 and 2 is too large for interactive computation") - # If not all genes, they used a var filter - if df1.var.shape[0] < self.gene_count: - mode = DiffExpMode.VAR_FILTER - if top_n: - raise Warning("Top N was specified but will not be used in 'Var Filter' mode") - else: - mode = DiffExpMode.TOP_N - if not top_n: - top_n = DEFAULT_TOP_N - - genes_idx = df1.var.index - # ensure we are using a dense ndarray - X1 = df1._X.toarray() if sparse.issparse(df1._X) else df1._X - X2 = df2._X.toarray() if sparse.issparse(df2._X) else df2._X - diffexp_result = stats.ttest_ind(X1, X2) - tstats = self._nan_to_zero(diffexp_result.statistic) - pval = self._nan_to_one(diffexp_result.pvalue) - bonferroni_pval = 1 - (1 - pval) ** self.gene_count - ave_exp_set1 = np.mean(X1, axis=0) - ave_exp_set2 = np.mean(X2, axis=0) - ave_diff = ave_exp_set1 - ave_exp_set2 - if mode == DiffExpMode.TOP_N: - sort_order = np.argsort(np.abs(tstats))[::-1] - # If top_n > length it will just return length - genes = self._top_sort(genes_idx, sort_order, top_n) - pval = self._top_sort(pval, sort_order, top_n) - bonferroni_pval = self._top_sort(bonferroni_pval, sort_order, top_n) - ave_exp_set1 = self._top_sort(ave_exp_set1, sort_order, top_n) - ave_exp_set2 = self._top_sort(ave_exp_set2, sort_order, top_n) - ave_diff = self._top_sort(ave_diff, sort_order, top_n) - - # varIndex, avgDiff, pVal, pValAdj, set1AvgExp, set2AvgExp - result = [] - for i in range(len(genes)): - result.append([genes[i], ave_diff[i], pval[i], bonferroni_pval[i], ave_exp_set1[i], ave_exp_set2[i]]) - # Results need to be returned in var index order - return sorted(result, key=lambda gene: gene[0]) + 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) + return sorted(result, key=lambda r: r[0]) def layout(self, filter, interactive_limit=None): """ @@ -404,8 +336,8 @@ class ScanpyEngine(CXGDriver): :return: [cellid, x, y, ...] """ try: - df = self.filter_dataframe(filter, include_uns=True) - except KeyError as e: + df = self.filter_dataframe(filter) + except (KeyError, IndexError) as e: raise FilterError(f"Error parsing filter: {e}") from e if interactive_limit and len(df.obs.index) > interactive_limit: raise InteractiveError("Size data is too large for interactive computation") diff --git a/server/test/test_api.py b/server/test/test_api.py index 9210ca73..7a0786c8 100644 --- a/server/test/test_api.py +++ b/server/test/test_api.py @@ -198,6 +198,7 @@ class EndPoints(unittest.TestCase): url = f"{URL_BASE}{endpoint}" params = { "mode": "topN", + "count": 10, "set1": { "filter": { "obs": { diff --git a/server/test/test_scanpy_engine.py b/server/test/test_scanpy_engine.py index 6861dd68..7ab0b971 100644 --- a/server/test/test_scanpy_engine.py +++ b/server/test/test_scanpy_engine.py @@ -85,7 +85,7 @@ class UtilTest(unittest.TestCase): } } } - data = self.data.filter_dataframe(filter_["filter"], include_uns=False) + data = self.data.filter_dataframe(filter_["filter"]) self.assertEqual(data.shape[1], 1) def test_filter_complex(self): @@ -184,7 +184,7 @@ class UtilTest(unittest.TestCase): layout = self.data.layout(filter_["filter"]) self.assertEqual(len(layout["coordinates"]), 497) - def test_diffexp(self): + def test_diffexp_topN(self): f1 = { "filter": { "obs": { @@ -199,11 +199,11 @@ class UtilTest(unittest.TestCase): } } } - result = self.data.diffexp(f1["filter"], f2["filter"]) + result = self.data.diffexp_topN(f1["filter"], f2["filter"]) self.assertEqual(len(result), 10) var_idx = [i[0] for i in result] self.assertEqual(var_idx, sorted(var_idx)) - result = self.data.diffexp(f1["filter"], f2["filter"], 20) + result = self.data.diffexp_topN(f1["filter"], f2["filter"], 20) self.assertEqual(len(result), 20) def test_data_frame(self):