gene sets summary route (#2099)

* gene sets summary route

* lint

* clarify return type

* style

local_server/app/app.py

@@ -155,6 +155,12 @@ class GenesetsAPI(Resource):
         return common_rest.genesets_put(request, data_adaptor)
 
 
+class GenesetSummaryAPI(Resource):
+    @rest_get_data_adaptor
+    def get(self, data_adaptor):
+        return common_rest.geneset_summary_get(request, data_adaptor)
+
+
 def get_api_base_resources(bp_base):
     """Add resources that are accessed from the api url"""
     api = Api(bp_base)
@@ -181,6 +187,7 @@ def get_api_dataroot_resources(bp_dataroot):
     add_resource(AnnotationsVarAPI, "/annotations/var")
     add_resource(DataVarAPI, "/data/var")
     add_resource(GenesetsAPI, "/genesets")
+    add_resource(GenesetSummaryAPI, "/geneset_summary")
     # Display routes
     add_resource(ColorsAPI, "/colors")
     # Computation routes

local_server/common/annotations/annotations.py

@@ -110,7 +110,7 @@ class Annotations(metaclass=ABCMeta):
         from io import StringIO
         import csv
 
-        if type(genesets) == dict:
+        if isinstance(genesets, dict):
             genesets = genesets.values()
 
         with StringIO() as sio:

local_server/common/annotations/local_file_csv.py

@@ -26,10 +26,14 @@ class AnnotationsLocalFile(Annotations):
         self.label_lock = threading.RLock()
         self.gene_sets_lock = threading.RLock()
 
-        # cache the most recent annotations.
+        # cache the most recent cell labels/annotations.
-        self.last_fname = None
+        self.last_label_fname = None
         self.last_labels = None
 
+        # cache the most recent gene sets.
+        self.last_geneset_fname = None
+        self.last_geneset = None
+
         # txn ID - used to de-dup geneset writes
         self.last_geneset_tid = 0
 
@@ -63,14 +67,14 @@ class AnnotationsLocalFile(Annotations):
         with self.label_lock:
             if fname is not None and os.path.exists(fname) and os.path.getsize(fname) > 0:
                 # returned the cached labels if possible, otherwise read them from the file
-                if fname == self.last_fname:
+                if fname == self.last_label_fname:
                     return self.last_labels
                 else:
                     labels = pd.read_csv(
                         fname, dtype="category", index_col=0, header=0, comment="#", keep_default_na=False
                     )
                     # update the cache
-                    self.last_fname = fname
+                    self.last_label_fname = fname
                     self.last_labels = labels
                     return labels
             else:
@@ -102,7 +106,7 @@ class AnnotationsLocalFile(Annotations):
                 open(fname, "w").close()
 
             # update the cache
-            self.last_fname = fname
+            self.last_label_fname = fname
             self.last_labels = df
 
     def read_gene_sets(self, data_adaptor, context=None):
@@ -116,8 +120,19 @@ class AnnotationsLocalFile(Annotations):
         with self.gene_sets_lock:
             tid = self.last_geneset_tid  # inside the critical section
             if fname is not None and os.path.exists(fname) and os.path.getsize(fname) > 0:
-                with open(fname, newline="") as f:
+                # return the cached genesets if possible, otherwise read from file and validate them
-                    gene_sets = read_gene_set_tidycsv(f, context)
+                if fname == self.last_geneset_fname:
+                    gene_sets = self.last_geneset
+                else:
+                    with open(fname, newline="") as f:
+                        gene_sets = read_gene_set_tidycsv(f, context)
+
+                    # validate
+                    gene_sets = data_adaptor.check_new_gene_sets(gene_sets, context)
+
+                    # update cache
+                    self.last_geneset_fname = fname
+                    self.last_geneset = gene_sets
 
         return (gene_sets, tid)
 
@@ -127,6 +142,9 @@ class AnnotationsLocalFile(Annotations):
         if type(tid) != int or tid < 0:
             raise ValueError("tid must be a positive integer")
 
+        # may raise
+        gene_sets = data_adaptor.check_new_gene_sets(gene_sets)
+
         with self.gene_sets_lock:
             # skip if the request is stale
             if tid is not None:
@@ -149,6 +167,10 @@ class AnnotationsLocalFile(Annotations):
                 f.write(header)
                 f.write(self.gene_sets_to_csv(gene_sets))
 
+            # update the cache
+            self.last_geneset_fname = fname
+            self.last_geneset = gene_sets if type(gene_sets) == dict else {g["geneset_name"]: g for g in gene_sets}
+
     def _get_userdata_idhash(self, data_adaptor):
         """
         Return a short hash that weakly identifies the user and dataset.

local_server/common/config/dataset_config.py

@@ -33,9 +33,9 @@ class DatasetConfig(BaseConfig):
                 "obo_location"
             ]
             self.user_annotations__gene_sets__readonly = default_config["user_annotations"]["gene_sets"]["readonly"]
-            self.user_annotations__local_file_csv__gene_sets_file = default_config["user_annotations"]["local_file_csv"][
+            self.user_annotations__local_file_csv__gene_sets_file = default_config["user_annotations"][
-                "gene_sets_file"
+                "local_file_csv"
-            ]
+            ]["gene_sets_file"]
 
             self.embeddings__names = default_config["embeddings"]["names"]
             self.embeddings__enable_reembedding = default_config["embeddings"]["enable_reembedding"]
@@ -169,10 +169,10 @@ class DatasetConfig(BaseConfig):
         if server_config.single_dataset__datapath:
             data_adaptor = self.get_data_adaptor()
             if self.user_annotations__local_file_csv__file:
-                data_adaptor.check_new_labels(self.user_annotations.read_labels(data_adaptor))
+                self.user_annotations.read_labels(data_adaptor)
             if self.user_annotations__local_file_csv__gene_sets_file:
                 try:
-                    data_adaptor.check_new_gene_sets(self.user_annotations.read_gene_sets(data_adaptor, context), context)
+                    self.user_annotations.read_gene_sets(data_adaptor, context)
                 except (ValueError, AnnotationsError, KeyError) as e:
                     raise ConfigurationError(f"Unable to read genesets CSV file: {str(e)}") from e
 

local_server/common/errors.py

@@ -55,4 +55,5 @@ define_exception("OntologyLoadFailure", "Raised when reading the ontology file f
 define_exception("ConfigurationError", "Raised when checking configuration errors")
 define_exception("PrepareError", "Raised when data is misprepared")
 define_exception("SecretKeyRetrievalError", "Raised when get_secret_key from AWS fails")
-define_exception("ObsoleteRequest", "Raised when the request is no longer valid.")
+define_exception("ObsoleteRequest", "Raised when the request is no longer valid.")
+define_exception("UnsupportedSummaryMethod", "Raised when a gene set summary method is unknown or unsupported.")

local_server/common/rest.py

@@ -19,6 +19,7 @@ from local_server.common.errors import (
     ColorFormatException,
     AnnotationsError,
     ObsoleteRequest,
+    UnsupportedSummaryMethod,
 )
 
 import json
@@ -336,7 +337,7 @@ def genesets_get(request, data_adaptor):
 
     try:
         annotations = data_adaptor.dataset_config.user_annotations
-        (genesets, tid) = data_adaptor.check_new_gene_sets(annotations.read_gene_sets(data_adaptor))
+        (genesets, tid) = annotations.read_gene_sets(data_adaptor)
 
         if preferred_mimetype == "text/csv":
             return make_response(
@@ -373,10 +374,39 @@ def genesets_put(request, data_adaptor):
         if genesets is None:
             abort(HTTPStatus.BAD_REQUEST)
 
-        (gs, _) = data_adaptor.check_new_gene_sets((genesets, tid))
+        annotations.write_gene_sets(genesets, tid, data_adaptor)
-        annotations.write_gene_sets(gs, tid, data_adaptor)
         return make_response(jsonify({"status": "OK"}), HTTPStatus.OK)
     except (ValueError, DisabledFeatureError, KeyError) as e:
         return abort_and_log(HTTPStatus.BAD_REQUEST, str(e), include_exc_info=True)
     except (ObsoleteRequest, TypeError) as e:
         return abort(HTTPStatus.NOT_FOUND, description=str(e))
+
+
+def geneset_summary_get(request, data_adaptor):
+    preferred_mimetype = request.accept_mimetypes.best_match(["application/octet-stream"])
+    if preferred_mimetype != "application/octet-stream":
+        return abort(HTTPStatus.NOT_ACCEPTABLE)
+
+    geneset_name = request.args.get("geneset_name", default=None)
+    summary_method = request.args.get("method", default="mean")
+    request_tid = request.args.get("tid", default=None)
+
+    try:
+        annotations = data_adaptor.dataset_config.user_annotations
+        (genesets, tid) = annotations.read_gene_sets(data_adaptor)
+
+        if request_tid is not None and int(request_tid) != tid:
+            return abort(HTTPStatus.NOT_FOUND, "Obsolete TID")
+        if geneset_name is None or geneset_name not in genesets:
+            return abort(HTTPStatus.BAD_REQUEST, "Gene set name not found.")
+        genes = [g["gene_symbol"] for g in genesets.get(geneset_name)["genes"]]
+
+        return make_response(
+            data_adaptor.get_gene_set_summary(geneset_name, genes, summary_method),
+            HTTPStatus.OK,
+            {"Content-Type": "application/octet-stream"},
+        )
+    except (ValueError) as e:
+        return abort(HTTPStatus.NOT_FOUND, description=str(e))
+    except (UnsupportedSummaryMethod) as e:
+        return abort(HTTPStatus.BAD_REQUEST, description=str(e))

local_server/data_anndata/anndata_adaptor.py

@@ -5,6 +5,7 @@ import anndata
 import numpy as np
 from packaging import version
 from pandas.core.dtypes.dtypes import CategoricalDtype
+import pandas as pd
 from scipy import sparse
 from server_timing import Timing as ServerTiming
 
@@ -12,7 +13,7 @@ import local_server.compute.diffexp_generic as diffexp_generic
 from local_server.common.colors import convert_anndata_category_colors_to_cxg_category_colors
 from local_server.common.constants import Axis, MAX_LAYOUTS
 from local_server.common.corpora import corpora_get_props_from_anndata
-from local_server.common.errors import PrepareError, DatasetAccessError, FilterError
+from local_server.common.errors import PrepareError, DatasetAccessError, FilterError, UnsupportedSummaryMethod
 from local_server.common.utils.type_conversion_utils import get_schema_type_hint_of_array
 from local_server.compute.scanpy import scanpy_umap
 from local_server.data_common.data_adaptor import DataAdaptor
@@ -367,3 +368,29 @@ class AnndataAdaptor(DataAdaptor):
     def get_var_keys(self):
         # return list of keys
         return self.data.var.keys().to_list()
+
+    def get_gene_set_summary(self, geneset_name, genes, method):
+        if method != "mean":
+            raise UnsupportedSummaryMethod("Unknown gene set summary method.")
+
+        var_index = self.parameters.get("var_names")
+        obs_selector, var_selector = self._filter_to_mask(
+            {
+                "var": {
+                    "annotation_value": [
+                        {
+                            "name": var_index,
+                            "values": genes,
+                        }
+                    ]
+                }
+            }
+        )
+
+        X = self.get_X_array(obs_selector, var_selector)
+        if sparse.issparse(X):
+            mean = X.mean(axis=1)
+        else:
+            mean = X.mean(axis=1, keepdims=True)
+        col_idx = pd.Index([geneset_name])
+        return encode_matrix_fbs(mean, col_idx=col_idx, row_idx=None)

local_server/data_common/data_adaptor.py

@@ -17,7 +17,7 @@ class DataAdaptor(metaclass=ABCMeta):
     """Base class for loading and accessing matrix data"""
 
     def __init__(self, data_locator, app_config, dataset_config=None):
-        if type(app_config) != AppConfig:
+        if not isinstance(app_config, AppConfig):
             raise TypeError("config expected to be of type AppConfig")
 
         # location to the dataset
@@ -157,7 +157,7 @@ class DataAdaptor(metaclass=ABCMeta):
     def _index_filter_to_mask(self, filter, count):
         mask = np.zeros((count,), dtype=np.bool)
         for i in filter:
-            if type(i) == list:
+            if isinstance(i, list):
                 mask[i[0] : i[1]] = True
             else:
                 mask[i] = True
@@ -262,15 +262,17 @@ class DataAdaptor(metaclass=ABCMeta):
 
         return labels_df
 
-    def check_new_gene_sets(self, args, context=None):
+    def check_new_gene_sets(self, genesets, context=None):
         """
         Check validity of gene sets, return if correct, else raise error.
         May also modify the gene set for conditions that should be resolved,
         but which do not warrant a hard error.
 
-        Argument 'args' must be a tuple containing (genesets, tid).  Genesets
+        Argument genesets may be either the REST OTA format (list of dicts) or the internal
-        may be either the REST OTA format (list of dicts) or the internal format
+        format (dict of dicts, keyed by the geneset name).
-        (dict of dicts, keyed by the geneset name).
+
+        Will return a modified genesets (eg, remove dups) of the same type as the
+        provided argument. Ie, dict->dict, list->list
 
         Rules:
         0. all geneset names must be unique.
@@ -283,17 +285,16 @@ class DataAdaptor(metaclass=ABCMeta):
            will generate a warning and the symbol removed.
         3. Duplicate gene symbols are silently de-duped.
         """
-        (genesets, tid) = args
         messagefn = context["messagefn"] if context else (lambda x: None)
 
         # accept genesets args as either the internal (dict) or REST (list) format,
         # as they are identical except for the dict being keyed by geneset_name.
-        if type(genesets) not in (dict, list):
+        if not isinstance(genesets, dict) and not isinstance(genesets, list):
             raise ValueError("Genesets must be either dict or list.")
-        genesets = genesets if type(genesets) == list else genesets.values()
+        genesets_iterable = genesets if isinstance(genesets, list) else genesets.values()
 
         # 0. check for uniqueness of geneset names
-        geneset_names = [gs["geneset_name"] for gs in genesets]
+        geneset_names = [gs["geneset_name"] for gs in genesets_iterable]
         if len(set(geneset_names)) != len(geneset_names):
             raise KeyError("All geneset names must be unique.")
 
@@ -316,18 +317,19 @@ class DataAdaptor(metaclass=ABCMeta):
         # 2. & 3. check for duplicate gene symbols, and those not present in the dataset. They will
         # generate a warning and be removed.
         var_names = set(self.query_var_array(self.parameters.get("var_names")))
-        for geneset in genesets:
+        for geneset in genesets_iterable:
-            if type(geneset) != dict:
+            if not isinstance(geneset, dict):
                 raise ValueError("Each geneset must be a dict.")
             geneset_name = geneset["geneset_name"]
             genes = geneset["genes"]
-            if type(genes) != list:
+            if not isinstance(genes, list):
                 raise ValueError("Geneset genes field must be a list")
+            geneset.setdefault("geneset_description", "")
             gene_symbol_already_seen = set()
             new_genes = []
             for gene in genes:
                 gene_symbol = gene["gene_symbol"]
-                if type(gene_symbol) != str or len(gene_symbol) == 0:
+                if not isinstance(gene_symbol, str) or len(gene_symbol) == 0:
                     raise ValueError("Gene symbol must be non-null string.")
                 if gene_symbol in gene_symbol_already_seen:
                     # duplicate check
@@ -345,11 +347,12 @@ class DataAdaptor(metaclass=ABCMeta):
                     continue
 
                 gene_symbol_already_seen.add(gene_symbol)
+                gene.setdefault("gene_description", "")
                 new_genes.append(gene)
 
             geneset["genes"] = new_genes
 
-        return args
+        return genesets
 
     def data_frame_to_fbs_matrix(self, filter, axis):
         """
@@ -479,3 +482,7 @@ class DataAdaptor(metaclass=ABCMeta):
         except RuntimeError:
             lastmod = None
         return lastmod
+
+    @abstractmethod
+    def get_gene_set_summary(self, geneset_name, genes, method):
+        pass

local_server/test/fixtures/pbmc3k-genesets.csv

@@ -10,3 +10,7 @@ second gene set,,SIK1
 third gene set,,NO_SUCH_GENE
 fourth_gene_set,fourth description,,gene intentionally missing
 fifth_dataset,,,
+summary test,,ACD,
+summary test,,AATF,
+summary test,,F5,
+summary test,,PIGU,

local_server/test/unit/common/test_api.py

@@ -509,6 +509,16 @@ class EndPointsAnnDataGenesets(unittest.TestCase, EndPoints):
                     {"genes": [], "geneset_description": "", "geneset_name": "third gene set"},
                     {"genes": [], "geneset_description": "fourth description", "geneset_name": "fourth_gene_set"},
                     {"genes": [], "geneset_description": "", "geneset_name": "fifth_dataset"},
+                    {
+                        "genes": [
+                            {"gene_description": "", "gene_symbol": "ACD"},
+                            {"gene_description": "", "gene_symbol": "AATF"},
+                            {"gene_description": "", "gene_symbol": "F5"},
+                            {"gene_description": "", "gene_symbol": "PIGU"},
+                        ],
+                        "geneset_description": "",
+                        "geneset_name": "summary test",
+                    },
                 ],
                 "tid": 0,
             },
@@ -531,6 +541,10 @@ second gene set,,SIK1,\r
 third gene set,,,\r
 fourth_gene_set,fourth description,,\r
 fifth_dataset,,,\r
+summary test,,ACD,\r
+summary test,,AATF,\r
+summary test,,F5,\r
+summary test,,PIGU,\r
 """,
         )
 
@@ -677,8 +691,86 @@ fifth_dataset,,,\r
             original_data,
         )
 
-    """
+    def test_get_geneset_summary(self):
-    TODO once we have some code to support it:
+        endpoint = "geneset_summary?geneset_name=summary%20test&method=mean"
-    1. GET genesets_summary
+        url = f"{self.URL_BASE}{endpoint}"
-    2. genesets_summary obeys tid
+        header = {"Accept": "application/octet-stream"}
-    """
+        result = self.session.get(url, headers=header)
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+        self.assertEqual(result.headers["Content-Type"], "application/octet-stream")
+        df = decode_fbs.decode_matrix_FBS(result.content)
+        self.assertEqual(df["n_rows"], 2638)
+        self.assertEqual(df["n_cols"], 1)
+        self.assertEqual(df["col_idx"], ["summary test"])
+        self.assertAlmostEqual(df["columns"][0][0], -0.19863907)
+
+    def test_get_geneset_summary_default_method(self):
+        endpoint = "geneset_summary?geneset_name=summary%20test"
+        url = f"{self.URL_BASE}{endpoint}"
+        header = {"Accept": "application/octet-stream"}
+        result = self.session.get(url, headers=header)
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+        self.assertEqual(result.headers["Content-Type"], "application/octet-stream")
+        df = decode_fbs.decode_matrix_FBS(result.content)
+        self.assertEqual(df["n_rows"], 2638)
+        self.assertEqual(df["n_cols"], 1)
+        self.assertEqual(df["col_idx"], ["summary test"])
+        self.assertAlmostEqual(df["columns"][0][0], -0.19863907)
+
+    def test_get_geneset_summary_check_tid(self):
+        # get the TID
+        result = self.session.get(f"{self.URL_BASE}genesets", headers={"Accept": "application/json"})
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+        tid = result.json()["tid"]
+
+        # current tid
+        endpoint = f"geneset_summary?geneset_name=summary%20test&tid={tid}"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+
+        # future tid
+        endpoint = f"geneset_summary?geneset_name=summary%20test&tid={tid+1}"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.NOT_FOUND)
+
+        # past tid
+        endpoint = f"geneset_summary?geneset_name=summary%20test&tid={tid-1}"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.NOT_FOUND)
+
+        # No tid - ie, skip check
+        endpoint = "geneset_summary?geneset_name=summary%20test"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+
+    def test_get_geneset_summary_edge_cases(self):
+        # attempt to summarize _all_ genesets, including edge cases with zero or one gene
+        result = self.session.get(f"{self.URL_BASE}genesets", headers={"Accept": "application/json"})
+        self.assertEqual(result.status_code, HTTPStatus.OK)
+        geneset_names = [gs["geneset_name"] for gs in result.json()["genesets"]]
+
+        for gs in geneset_names:
+            endpoint = f"geneset_summary?geneset_name={gs}"
+            result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+            self.assertEqual(result.status_code, HTTPStatus.OK)
+            self.assertEqual(result.headers["Content-Type"], "application/octet-stream")
+            df = decode_fbs.decode_matrix_FBS(result.content)
+            self.assertEqual(df["n_rows"], 2638)
+            self.assertEqual(df["n_cols"], 1)
+            self.assertEqual(df["col_idx"], [gs])
+
+    def test_get_geneset_error_handling(self):
+        # no geneset
+        endpoint = "geneset_summary"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.BAD_REQUEST)
+
+        # unknown geneset
+        endpoint = "geneset_summary?geneset_name=NO_SUCH_GENE_SET"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.BAD_REQUEST)
+
+        # unknown method
+        endpoint = "geneset_summary?geneset_name=summary%20test&method=NO_SUCH_METHOD"
+        result = self.session.get(f"{self.URL_BASE}{endpoint}", headers={"Accept": "application/octet-stream"})
+        self.assertEqual(result.status_code, HTTPStatus.BAD_REQUEST)