diff --git a/MANIFEST.in b/MANIFEST.in index 7d059a83..90913b0c 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -3,3 +3,5 @@ recursive-include server/common/web/static * include server/requirements.txt include server/requirements-prepare.txt +include server/converters/schema/hgnc_complete_set.txt.gz +include server/converters/schema/schema_definitions * diff --git a/dev_docs/schema_guide.md b/dev_docs/schema_guide.md new file mode 100644 index 00000000..cd70611a --- /dev/null +++ b/dev_docs/schema_guide.md @@ -0,0 +1,175 @@ +# Cellxgene Schema Guide + +Datasets included in the [data portal](https://cellxgene.cziscience.com/) and hosted cellxgene need to follow the schema +described [here](https://github.com/chanzuckerberg/corpora-data-portal/blob/main/backend/schema/corpora_schema.md). That +schema defines some required fields, requirements about feature labels, and some optional fields that mostly help with +presentation. + +The number of fields is rather low, and we expect that information needed to populate those fields should either already +be present in datasets prepared by a submitter or be easy to obtain. However, this still leaves the task of actually +manipulating the dataset so that it follows the schema: adjusting field names, ensuring proper ontologies are used, +converting gene symbols to a common set, etc. This can be tedious and error-prone, and at the beginning of the hosted +cellxgene project, this was always done with engineering support. As we increase the rate at which we add data, we want +to eliminate the need for engineering support so that ultimately submitters themselves can create files that follow the +schema. + +## `cellxgene schema apply` + +To enable this, we have a new cellxgene subcommand, `cellxgene schema`, that handles applying and verifying the schema. +Its first subcommand, `cellxgene schema apply`, takes three inputs: + +1. A source h5ad file. The input needs to be an AnnData file, so if a submitter has, say, a serialized Seurat or + SingleCellExperiment object, it needs to be converted to AnnData first. This can be done with + [sceasy](https://github.com/cellgeni/sceasy) or via + [Seurat](https://satijalab.org/seurat/v3.1/conversion_vignette.html). +2. A configuration yaml file that describes the fields to add and conversions to apply (see below). +3. A name for the new h5ad file that should follow the schema. + +### Configuration yaml + +The configuration yaml file describes how to apply the schema. This is an example of a "skeleton" yaml that has all the +fields required for the 1.0.0 schema but is not yet filled in with any logic: + +``` +uns: + version: + corpora_schema_version: 1.0.0 + corpora_encoding_version: 0.1.0 + contributors: + title: + layer_descriptions: + preprint_doi: + publication_doi: + organism_ontology_term_id: +obs: + tissue_ontology_term_id: + assay_ontology_term_id: + disease_ontology_term_id: + cell_type_ontology_term_id: + sex: + ethnicity_ontology_term_id: + development_stage_ontology_term_id: +fixup_gene_symbols: +``` + +#### Unstructured metadata +The first section is `uns`, which includes metadata fields that describe the whole dataset (see +[here](https://anndata.readthedocs.io/en/latest/) for further description of `uns` and `obs`.). + +The first line is `version`, which is required for most of our tooling to work. The schema version is set at +1.0.0 in the example above, but of course for future versions that should be changed. + +Next is `contributors` which describes who is adding the dataset to the portal. If you consult the schema, you see that +contributors is a list where each element can have `name`, `email`, and `institution`. So when filled out, the +`contributors` field should look like this: + +``` +contributors: + - name: Mary B. Scientist + email: mbs@singlecell.edu + institution: Single-Cell University + - name: Robert J. Scientist + email: rjs@usingle.edu + institution: University of Single Cell +``` + +`title` is the name of the dataset, and is just a string that gets displayed in the portal and cellxgene to identify the +dataset. + +`layer_descriptions` is free text descriptions of the different +[layers](https://anndata.readthedocs.io/en/latest/anndata.AnnData.layers.html) of the AnnData file. It should look like +this when complete, depending on what layers are present: +``` +layer_descriptions: + X: CPM and logged + raw.X: raw +``` +Note that one of the layers needs to be "raw", that is, the AnnData file must contain raw counts. + +The two DOI fields are optional but can be included if the dataset is associated with a publication or preprint. Note +that the DOI should be a full url: +``` +publication_doi: https://doi.org/10.1073%2Fpnas.83.15.5372 +``` + +Finally, the `organism_ontology_term_id` field is the species of the donor organism from the NCBITaxon ontology. The +value for _Homo sapiens_ is `NCBITaxon:9606`: +``` +organism_ontology_term_id: NCBITaxon:9606 +``` +Note that the schema also requires a human-readable `organism` field, but this doesn't need to be included in the yaml. +When the `cellxgene schema apply` script encounters an ontology field, it looks up the label for the term(s) and inserts it +into the appropriate field. + + +#### Observation metadata +The next section is `obs`, which is metadata than can vary for each observation (and "observation" usually means cell). +These fields are all ontology fields except for `sex`, which has its own enumerated set of permitted values. + +There are two ways to fill in the `obs` fields. The first is useful when there is only one value for all the +observations in the dataset. This is not uncommon, for example all cells often come from the same assay. In that case +just insert the ontology term: +``` +assay_ontology_term_id: EFO:0009922 +``` + +The second is for when there is an existing field in the dataset that needs to be mapped to the schema field. For +example, the submitter may have included cell type annotations in a field called `CellType`, and those annotations may +just be free text. This doesn't follow the schema because it needs to be in `cell_type_ontology_term_id` and +`cell_type`, and it needs ontology terms and labels, not just any text. In that case the field can be a dictionary: + +``` +cell_type_ontology_term_id: + CellType: + t-cell: CL:0000084 + b-cell: CL:0000236 +``` + +This will look at the `obs.CellType` field in the dataset, and where it has the value "t-cell", it will insert +`CL:0000084` into `cell_type_ontology_term_id` and its label `T cell` into `cell_type`. + +Now there are often situations where there is no valid ontology term for some field. For example, the dataset may have +been produced via an assay not present in `EFO`. Or, a particular cell type may have no entry in `CL`. In that case, a +free text description can be used in the `ontology_term_id` field: + +``` +assay_ontology_term_id: Sci-Plex +cell_type_ontology_term_id: + CellType: + t-cell: CL:0000084 + b-cell: CL:0000236 + new cell type: new cell type +``` + +In these cases, the `cellxgene schema apply` script will leave the ontology field blank and move the free text +description into the label field. So the `assay_ontology_term_id` in the new dataset would be `""` but `assay` would be +`Sci-Plex`. + + +#### Gene symbol harmonization + +The last section describes how gene symbol conversion should be applied to each of the layers. This is similar to the +`layer_descriptions` field above, but there are only three permitted values: `raw`, `log1p`, and `sqrt`: + +``` +fixup_gene_symbols: + X: log1p + raw.X: raw +``` + +This tells the script how each each layer was transformed from raw values that can be directly summed. `raw` means that +the layer contains raw counts or some linear tranformation of raw counts. `log1p` means that the layer has `log(X + 1)` +for each the raw `X` values. `sqrt` means `sqrt(X)` (this is not common). For layers produced by Seurat's normalization +or SCTransform functions, the correct choice is usually `log1p`. + + +### `cellxgene schema validate` + +The next `cellxgene schema` subcommand is `cellxgene schema validate`, and it validates that a given h5ad follows a +version of the schema. It accepts two parameters: + +1. The h5ad file to check +2. The version of the schema to check against. + +If the validation succeeds, the command will have a zero exit code. If it does not, it will have a non-zero exit code +and will print validation failure messages. diff --git a/server/cli/cli.py b/server/cli/cli.py index dc3e5837..f8f7fde6 100644 --- a/server/cli/cli.py +++ b/server/cli/cli.py @@ -4,6 +4,7 @@ from .convert_to_cxg import convert_to_cxg from .launch import launch from .prepare import prepare from .upgrade import log_upgrade_check +from .schema import schema_cli from .. import __version__ @@ -31,3 +32,4 @@ def cli(upgrade_check): cli.add_command(launch) cli.add_command(prepare) cli.add_command(convert_to_cxg) +cli.add_command(schema_cli) diff --git a/server/cli/schema.py b/server/cli/schema.py new file mode 100644 index 00000000..5f16ec64 --- /dev/null +++ b/server/cli/schema.py @@ -0,0 +1,72 @@ +import click + +from server.converters.schema import remix, validate + + +@click.group( + name="schema", + subcommand_metavar="COMMAND ", + short_help="Apply and validate the cellxgene data integration schema to an h5ad file.", + context_settings=dict(max_content_width=85, help_option_names=["-h", "--help"]), +) +def schema_cli(): + try: + import scanpy # noqa: F401 + except ImportError: + raise click.ClickException( + "[cellxgene] cellxgene schema requires scanpy" + ) + + +@click.command( + name="apply", + short_help="(experimental) Apply the cellxgene data integration schema to an h5ad.", + help="(experimental) Using a yaml file that describes schema values to insert or convert and in input " + "h5ad file, apply the schema changes and create a new, conforming h5ad.", +) +@click.option( + "--source-h5ad", + help="Input h5ad file.", + nargs=1, + required=True, + type=click.Path(exists=True, dir_okay=False), +) +@click.option( + "--remix-config", + help="Config yaml with information on how to apply the schema.", + nargs=1, + required=True, + type=click.Path(exists=True, dir_okay=False), +) +@click.option( + "--output-filename", + help="Filename for the new, schema-conforming h5ad file.", + required=True, + nargs=1 +) +def schema_apply(source_h5ad, remix_config, output_filename): + remix.apply_schema(source_h5ad, remix_config, output_filename) + + +@click.command( + name="validate", + short_help="(experimental) Check that an h5ad follows the cellxgene data integration schema.", +) +@click.argument( + "h5ad", + nargs=1, + type=click.Path(exists=True, dir_okay=False), +) +@click.option( + "--shallow", + help="When true, just check that the correct version information is present.", + default=False, + show_default=True, + is_flag=True, +) +def schema_validate(h5ad, shallow): + validate.validate(h5ad, shallow) + + +schema_cli.add_command(schema_apply) +schema_cli.add_command(schema_validate) diff --git a/server/converters/schema/__init__.py b/server/converters/schema/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/server/converters/schema/gene_symbol.py b/server/converters/schema/gene_symbol.py new file mode 100644 index 00000000..00dd386c --- /dev/null +++ b/server/converters/schema/gene_symbol.py @@ -0,0 +1,208 @@ +"""Helpers for converting and checking HGNC gene symbols.""" + +import argparse +import enum +import logging +import os +import re +import numpy as np +import pandas as pd + + +def get_upgraded_var_index(var, hgnc_path=None): + """Given an anndata var dataframe, return a new index for the dataframe + where human gene symbols have been upgraded to the current HGNC set. + """ + + if not hgnc_path: + hgnc_path = os.path.join(os.path.dirname(os.path.realpath(__file__)), "hgnc_complete_set.txt.gz") + + hgnc_symbol_checker = HGNCSymbolChecker.from_hgnc_records(hgnc_path) + + return pd.Index([hgnc_symbol_checker.upgrade_symbol(s) for s in var.index]) + + +class SymbolStatus(enum.Enum): + """The status of a symbol in the HGNC database. + + APPROVED: Currently a valid symbol + WITHDRAWN: A previously approved HGNC symbol for a gene that has since been shown + not to exist _unless_ that symbol is also approved + AMBIGUOUS: A symbol that is not approved but is an alias or previous symbol for + multiple approved symbols + UPGRADABLE: A symbol that is not approved but unambiguously maps to an approved + symbol + UNKNOWN: A symbol that does not appear in HGNC + """ + + APPROVED = 1 + WITHDRAWN = 2 + AMBIGUOUS = 3 + UPGRADABLE = 4 + UNKNOWN = 5 + + +class HGNCSymbolChecker: + """Handle checking and correcting HGNC symbols.""" + + def __init__(self, approved_symbols, withdrawn_symbols, ambiguous_symbols, symbol_map): + self.approved_symbols = approved_symbols + self.withdrawn_symbols = withdrawn_symbols + self.ambiguous_symbols = ambiguous_symbols + self.symbol_map = symbol_map + + def print_symbol_map(self): + """Print out a map from old symbol to new symbol.""" + + for symbol_pair in self.symbol_map.items(): + print("\t".join(symbol_pair)) + + def check_symbol(self, symbol): + """See if a symbol if approved or something else.""" + if symbol in self.approved_symbols: + return SymbolStatus.APPROVED + + if symbol in self.withdrawn_symbols: + return SymbolStatus.WITHDRAWN + + if symbol in self.ambiguous_symbols: + return SymbolStatus.AMBIGUOUS + + if symbol in self.symbol_map: + return SymbolStatus.UPGRADABLE + + return SymbolStatus.UNKNOWN + + def upgrade_symbol(self, symbol): + """Return the approved symbol for the given symbol. + + If the symbol cannot be upgraded, just return the original symbol. + """ + + fixed_symbol, stripped_symbol = format_symbol(symbol) + + if fixed_symbol in self.approved_symbols: + return fixed_symbol + elif fixed_symbol in self.symbol_map: + return self.symbol_map[fixed_symbol] + elif stripped_symbol in self.approved_symbols: + return stripped_symbol + elif stripped_symbol in self.symbol_map: + return self.symbol_map[stripped_symbol] + + return symbol + + @classmethod + def from_hgnc_records(cls, hgnc_dataset_path): + """Parse a hgnc database download into a HGNCSymbolChecker object.""" + + def all_symbols(record): + """Get all the symbols associated with an HGNC record including previous, alias, + and approved.""" + yield format_symbol(record["symbol"])[0] + for symbol in alias_and_previous_symbols(record): + yield symbol + + def alias_and_previous_symbols(record): + """Get alias and previous symbols from an HGNC record.""" + for field in ("alias_symbol", "prev_symbol"): + if record[field] is not np.nan: + for symbol in record[field].split("|"): + yield format_symbol(symbol)[0] + + hgnc_records = pd.read_csv(hgnc_dataset_path, sep="\t", header=0, low_memory=False).to_dict("records") + + # Get all symbols that are currently approved. + approved_symbols = set() + for record in hgnc_records: + if record["status"] == "Approved": + approved_symbols.add(format_symbol(record["symbol"])[0]) + + # Get all symbols that have been withdrawn + withdrawn_symbols = set() + for record in hgnc_records: + if record["status"] == "Entry Withdrawn": + for symbol in all_symbols(record): + withdrawn_symbols.add(symbol) + + # If a symbol is both approved and withdrawn, be optimistic and call it approved + logging.warning( + f"Some symbols are simulaneously withdrawn and approved\n" + f"We will treat them at approved:\n" + f"{withdrawn_symbols.intersection(approved_symbols)}" + ) + withdrawn_symbols = withdrawn_symbols.difference(approved_symbols) + + # Now try to map from symbols that are not approved but are an alias or previous symbol for an approved symbol + alias_previous_to_approved = {} + ambiguous_symbols = set() + + for record in hgnc_records: + if record["status"] == "Approved": + + # The approved symbol is what we'll map to + approved_symbol = format_symbol(record["symbol"])[0] + + for symbol in alias_and_previous_symbols(record): + + # If the alias or previous symbol is also an approved symbol, + # we'll just leave it alone + if symbol in approved_symbols: + continue + + # If the alias or previous symbol maps to a different approved symbol, mark it as ambiguous + if symbol in alias_previous_to_approved and alias_previous_to_approved[symbol] != approved_symbol: + ambiguous_symbols.add(symbol) + else: + alias_previous_to_approved[symbol] = approved_symbol + + # Remove all the ambiguous symbols from the map + for ambiguous_symbol in ambiguous_symbols: + alias_previous_to_approved.pop(ambiguous_symbol) + + return HGNCSymbolChecker(approved_symbols, withdrawn_symbols, ambiguous_symbols, alias_previous_to_approved) + + +def format_symbol(symbol): + """HGNC rules say symbols should all be upper case except for C#orf#. However, case is + variable in both alias and previous symbols as well as in the symbols we get in + submissions. So, upper case everything except for the one situation where mixed-case + is allowed, which are the genes like C2orf157. + + Also, seurat and scanpy append ".1" or "-1" to duplicated gene names, and these altered + names persist throughout the life of the object. They won't match against the HGNC database + and we want to merge them, so we need to strip off the suffix and try matching again. + + This function takes a symbol and returns the symbol with the fixed case and also with the + seurat/scanpy suffix stripped off. + """ + + match = re.match(r"^(C)(\d+)(orf)(\d+)$", symbol, re.IGNORECASE) + + if match: + fixed_case = f"C{match.group(2)}orf{match.group(4)}" + else: + fixed_case = symbol.upper() + + suffix_stripped = re.sub(r"[\.\-]\d+$", "", fixed_case) + + return fixed_case, suffix_stripped + + +def main(): + """When called as main, parse a given hgnc download and print out a map from old to new + symbol. + """ + parser = argparse.ArgumentParser() + parser.add_argument( + "hgnc_dataset", help="HGNC dataset tsv, available from www.genenames.org/download/statistics-and-files/" + ) + args = parser.parse_args() + + hgnc_symbol_checker = HGNCSymbolChecker.from_hgnc_records(args.hgnc_dataset) + + hgnc_symbol_checker.print_symbol_map() + + +if __name__ == "__main__": + main() diff --git a/server/converters/schema/hgnc_complete_set.txt.gz b/server/converters/schema/hgnc_complete_set.txt.gz new file mode 100644 index 00000000..29c3c7a9 Binary files /dev/null and b/server/converters/schema/hgnc_complete_set.txt.gz differ diff --git a/server/converters/schema/ontology.py b/server/converters/schema/ontology.py new file mode 100644 index 00000000..8a524402 --- /dev/null +++ b/server/converters/schema/ontology.py @@ -0,0 +1,86 @@ +"""Methods for working with ontologies and the OLS.""" +from urllib.parse import quote_plus + +import requests + +OLS_API_ROOT = "http://www.ebi.ac.uk/ols/api" + +# Curie means something like CL:0000001 + + +def _ontology_name(curie): + """Get the name of the ontology from the curie, CL or UBERON for example.""" + return curie.split(":")[0] + + +def _ontology_value(curie): + """Get the id component of the curie, 0000001 from CL:0000001 for example.""" + return curie.split(":")[1] + + +def _double_encode(url): + """Double url encode a url. This is required by the OLS API.""" + return quote_plus(quote_plus(url)) + + +def _iri(curie): + """Get the iri from a curie. This is a bit hopeful that they all map to purl.obolibrary.org""" + if _ontology_name(curie) == "EFO": + return f"http://www.ebi.ac.uk/efo/EFO_{_ontology_value(curie)}" + return f"http://purl.obolibrary.org/obo/{_ontology_name(curie)}_{_ontology_value(curie)}" + + +class OntologyLookupError(Exception): + """Exception for some problem with looking up ontology information.""" + + +def _ontology_info_url(curie): + """Get the to make a GET to to get information about an ontology term.""" + + # If the curie is empty, just return an empty string. This happens when there is no + # valid ontology value. + if not curie: + return "" + else: + return f"{OLS_API_ROOT}/ontologies/{_ontology_name(curie)}/terms/{_double_encode(_iri(curie))}" + + +def get_ontology_label(curie): + """For a given curie like 'CL:1000413', get the label like 'endothelial cell of artery'""" + + url = _ontology_info_url(curie) + + if not url: + return "" + + response = requests.get(url) + + if not response.ok: + raise OntologyLookupError( + f"Curie {curie} lookup failed, got status code {response.status_code}: {response.text}" + ) + return response.json()["label"] + + +def lookup_candidate_term(label, ontology="cl", method="select"): + """Lookup candidate terms for a label. This is useful when there is an existing label in a + submitted dataset, and you want to find an appropriate ontology term. + + Args: + label: the label to find ontology terms for + ontology: the ontology to search in, cl or uberon or efo for example + method: select or search. search provides much broader results + + Returns: + list of (curie, label) tuples returned by OLS + """ + # using OLS REST API [https://www.ebi.ac.uk/ols/docs/api] + url = f"{OLS_API_ROOT}/{method}?q={quote_plus(label)}&ontology={ontology.lower()}" + response = requests.get(url) + + if not response.ok: + raise OntologyLookupError( + f"Label {label} lookup failed, got status code {response.status_code}: {response.text}" + ) + + return [(r["obo_id"], r["label"]) for r in response.json()["response"]["docs"]] diff --git a/server/converters/schema/remix.py b/server/converters/schema/remix.py new file mode 100644 index 00000000..5cc60746 --- /dev/null +++ b/server/converters/schema/remix.py @@ -0,0 +1,264 @@ +import argparse +import collections +import json +import logging +import math +import string + +import anndata +import numpy as np +import pandas as pd +import yaml + +from . import gene_symbol +from . import ontology +from . import validate + +REPLACE_SUFFIX = "_original" +ONTOLOGY_SUFFIX = "_ontology_term_id" + + +def is_curie(value): + """Return True iff the value is an OBO-id CURIE like EFO:000001""" + return (value.count(":") + and all(len(part) > 0 for part in value.split(":")) + and all(c in string.digits for c in value.split(":")[1])) + + +def is_ontology_field(field_name): + """Return True iff the field_name is an ontology field like tissue_ontology_term_id""" + return field_name.endswith(ONTOLOGY_SUFFIX) + + +def get_label_field_name(field_name): + """Get the associated label field from an ontology field, assay_ontology_term_id --> assay""" + return field_name[: -len(ONTOLOGY_SUFFIX)] + + +def split_suffix(maybe_curie): + """Split off the (cell culture) or (organoid) suffix.""" + + suffixes = [" (cell culture)", " (organoid)"] + for suffix in suffixes: + if maybe_curie.endswith(suffix): + return maybe_curie[:-len(suffix)], suffix + return maybe_curie, "" + + +def get_curie_and_label(maybe_curie): + """Given a string that might be a curie, return a (curie, label) pair""" + + maybe_curie, suffix = split_suffix(maybe_curie) + if not is_curie(maybe_curie): + return ("", maybe_curie + suffix) + return (maybe_curie + suffix, ontology.get_ontology_label(maybe_curie) + suffix) + + +def safe_add_field(adata_attr, field_name, field_value): + """Add a field and value to an AnnData, but don't clobber an exising value.""" + + if ( + isinstance(field_value, list) + and field_value + and isinstance(field_value[0], dict) + ): + field_value = json.dumps(field_value) + if field_name in adata_attr: + adata_attr[field_name + REPLACE_SUFFIX] = adata_attr[field_name] + adata_attr[field_name] = field_value + + +def remix_uns(adata, uns_config): + """Add fields from the config to adata.uns""" + for field_name, field_value in uns_config.items(): + + if is_ontology_field(field_name): + # If it's an ontology field, look it up + label_field_name = get_label_field_name(field_name) + ontology_term, ontology_label = get_curie_and_label(field_value) + safe_add_field(adata.uns, field_name, ontology_term) + safe_add_field(adata.uns, label_field_name, ontology_label) + else: + safe_add_field(adata.uns, field_name, field_value) + + +def remix_obs(adata, obs_config): + """Add fields from the config to adata.obs""" + + for field_name, field_value in obs_config.items(): + + if isinstance(field_value, dict): + # If the value is a dict, that means we are supposed to map from an + # existing column to the new one + source_column, column_map = next(iter(field_value.items())) + nan_value = None + for key in column_map: + if isinstance(key, float) and math.isnan(key): + nan_value = column_map[key] + if nan_value is not None: + column_map["nan"] = nan_value + + for key in column_map: + if key not in adata.obs[source_column].unique(): + logging.warning(f'Key {key} not in adata.obs["{source_column}"]') + + for value in adata.obs[source_column].unique(): + if value not in column_map: + logging.warning(f'Value {value} in adata.obs["{source_column}"] not in translation dict') + + if is_ontology_field(field_name): + ontology_term_map, ontology_label_map = {}, {} + logging.info(f"Looking up labels for {field_name}") + for original_value, maybe_curie in column_map.items(): + curie, label = get_curie_and_label(maybe_curie) + ontology_term_map[original_value] = curie + ontology_label_map[original_value] = label + logging.info(f"Mapping {original_value} -> {curie} -> {label}") + + ontology_column = adata.obs[source_column].replace( + ontology_term_map, inplace=False + ) + label_column = adata.obs[source_column].replace( + ontology_label_map, inplace=False + ) + + safe_add_field(adata.obs, field_name, ontology_column) + safe_add_field( + adata.obs, get_label_field_name(field_name), label_column + ) + else: + label_column = adata.obs[source_column].replace( + column_map, inplace=False + ) + safe_add_field(adata.obs, field_name, label_column) + + else: + if is_ontology_field(field_name): + # If it's an ontology field, look it up + label_field_name = get_label_field_name(field_name) + ontology_term, ontology_label = get_curie_and_label(field_value) + safe_add_field(adata.obs, field_name, ontology_term) + safe_add_field(adata.obs, label_field_name, ontology_label) + else: + safe_add_field(adata.obs, field_name, field_value) + + +def merge_df(df, domain, index, columns): + """ + Given a dataframe with duplicate column labels, merge and return a dataframe where + the duplicates have been merged together, resulting in a dataframe with unique column + labels. + + "merge" depends on the value of domain. If the domain is "raw", then duplicate columns + can just be summed. If it's "log1p" or "sqrt", it needs to be exp1m'd or squared, then + summed, and then logged or sqrt'd again. + """ + + if not isinstance(df, np.ndarray): + to_merge = df.toarray() + else: + to_merge = df + if domain == "raw": + merged_df = pd.DataFrame(to_merge, index=index, columns=columns).sum( + axis=1, level=0, skipna=False + ) + elif domain == "log1p": + merged_df = ( + pd.DataFrame(np.expm1(to_merge, dtype=np.float128), index=index, columns=columns) + .sum(axis=1, level=0, skipna=False) + ) + merged_df = pd.DataFrame(np.log1p(merged_df.to_numpy()), index=merged_df.index, columns=merged_df.columns) + elif domain == "sqrt": + merged_df = ( + pd.DataFrame(np.square(to_merge), index=index, columns=columns) + .sum(axis=1, level=0, skipna=False) + ) + merged_df = pd.DataFrame(np.sqrt(merged_df.to_numpy()), index=merged_df.index, columns=merged_df.columns) + + return merged_df + + +def fixup_gene_symbols(adata, fixup_config): + """Update the var index to hold a consistent set of HGNC gene symbols.""" + + upgraded_var_index = gene_symbol.get_upgraded_var_index(adata.var) + + merged_X = merge_df(adata.X, fixup_config["X"], adata.obs.index, upgraded_var_index) + fixup_adata = anndata.AnnData( + X=merged_X, + obs=adata.obs, + var=merged_X.columns.to_frame(name="hgnc_gene_symbol"), + uns=adata.uns, + obsm=adata.obsm, + ) + + for layer, domain in fixup_config.items(): + if layer == "X": + continue + if layer == "raw.X": + df = adata.raw.X + else: + df = adata.layers[layer] + + merged_df = merge_df(df, domain, adata.obs.index, upgraded_var_index) + assert merged_df.index.equals(merged_X.index) + assert merged_df.columns.equals(merged_X.columns) + + if domain == "raw": + fixup_raw = anndata.AnnData( + X=merged_df, + obs=adata.obs, + var=merged_X.columns.to_frame(name="hgnc_gene_symbol"), + ) + fixup_adata.raw = fixup_raw + else: + fixup_adata.layers[layer] = merged_df + + return fixup_adata + +def _strip_version(adata): + """Remove version information from the AnnData object.""" + + if "version" in adata.uns_keys(): + del adata.uns["version"] + +def apply_schema(source_h5ad, remix_config, output_filename): + + try: + import scanpy + except ImportError: + raise ImportError("scanpy must be installed for cellxgene schema") + adata = scanpy.read_h5ad(source_h5ad) + config = yaml.load(open(remix_config), Loader=yaml.FullLoader) + remix_uns(adata, config["uns"]) + remix_obs(adata, config["obs"]) + + if config.get("fixup_gene_symbols"): + adata = fixup_gene_symbols(adata, config["fixup_gene_symbols"]) + + if ("version" in adata.uns_keys() + and isinstance(adata.uns["version"], collections.Mapping) + and "corpora_schema_version" in adata.uns["version"]): + schema_version = adata.uns["version"]["corpora_schema_version"] + try: + validate.get_schema_definition(schema_version) + except ValueError: + logging.warning(f"Stripping version information out of AnnData because schema " + f"version {schema_version} is unknown.") + _strip_version(adata) + + if not validate.validate_adata(adata, shallow=False): + logging.warning(f"Stripping version information out of AnnData because it does not " + f"follow schema version {schema_version} .") + _strip_version(adata) + + adata.write_h5ad(output_filename, compression="gzip") + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("--source-h5ad", required=True) + parser.add_argument("--remix-config", required=True) + parser.add_argument("--output-filename", required=True) + args = parser.parse_args() + apply_schema(args.source_h5ad, args.remix_config, args.output_filename) diff --git a/server/converters/schema/schema_definitions/1_0_0.yaml b/server/converters/schema/schema_definitions/1_0_0.yaml new file mode 100644 index 00000000..acff2238 --- /dev/null +++ b/server/converters/schema/schema_definitions/1_0_0.yaml @@ -0,0 +1,95 @@ +title: Corpora schema version 1.0.0 +type: anndata +components: + uns: + type: dict + keys: + version: + type: dict + keys: + corpora_schema_version: null + corpora_encoding_version: null + title: + type: string + contributors: + type: stringified list of dicts + layer_descriptions: + type: dict + keys: + X: null + organism: + type: string + nullable: false + organism_ontology_term_id: + type: curie + prefixes: + - NCBITaxon + var: + type: dataframe + index: + type: human-readable string + unique: true + obs: + type: dataframe + index: + unique: true + columns: + tissue: + type: human-readable string + nullable: false + tissue_ontology_term_id: + type: suffixed curie + nullable: true + prefixes: + - UBERON + assay: + type: human-readable string + nullable: false + assay_ontology_term_id: + type: curie + nullable: true + prefixes: + - EFO + disease: + type: human-readable string + nullable: false + disease_ontology_term_id: + type: curie + nullable: true + prefixes: + - MONDO + - PATO + cell_type: + type: human-readable string + nullable: false + cell_type_ontology_term_id: + type: curie + nullable: true + prefixes: + - CL + - UBERON + sex: + type: string + enum: + - male + - female + - mixed + - unknown + - other + ethnicity: + type: human-readable string + nullable: false + ethnicity_ontology_term_id: + type: curie + nullable: true + prefixes: + - HANCESTRO + development_stage: + type: human-readable string + nullable: false + development_stage_ontology_term_id: + type: curie + nullable: true + prefixes: + - HsapDv + - EFO diff --git a/server/converters/schema/validate.py b/server/converters/schema/validate.py new file mode 100644 index 00000000..ec463dd5 --- /dev/null +++ b/server/converters/schema/validate.py @@ -0,0 +1,236 @@ +import json +import re +import os +import sys + +import pandas as pd +import yaml + + +def _is_null(v): + """Return True if v is null, for one of the multiple ways a "null" value shows up in an h5ad.""" + return pd.isnull(v) or (hasattr(v, "__len__") and len(v) == 0) + + +def _validate_stringified_list_of_dicts(s): + """Verify that a string can be parsed into a list. + + We have some types that are lists of dicts. Those cannot be stored directly in an h5ad, so we have to + json.dumps them. This verifies that we can load them back. + """ + + try: + list_ = json.loads(s) + if not isinstance(list_, list): + return False + for el in list_: + if not isinstance(el, dict): + return False + return True + except (json.JSONDecodeError, TypeError): + pass + return False + + +def _validate_human_readable_string(s): + """Verify that a string is human-readable. + + There are parts of the schema where a "human-readable" string is required. "Human-readable" is kind + of vague and subjective. I feel like I can read many strings. So here we just check for the main ways + that fails: someone puts in an ontology term id or and ensembl gene/transcript id. + + Returns False if s is not a string or is one of those bad string types. + """ + + return isinstance(s, str) and (not re.match(r"[A-Z]\w+:\d+", s)) and (not re.match(r"ENS[GT]\d+$", s)) + + +def _validate_curie(c, prefixes): + """Verify that a string is a valid compact URI, like EFO:000001. If prefixes is not empty, make sure the + prefix of the curies is in prefixes. + """ + + if not c: + return True + + match = re.match(r"([A-Z]\w+):\d+$", c) + + if prefixes: + return match and match.group(1) in prefixes + else: + return match + + +def _validate_suffixed_curie(c, prefixes): + """Verify that a string is a compact URI with an optional suffix like 'EFO:00001 (cell culture)'""" + + # Pull off the suffix + suffix = re.findall(r"\ \(.*\)$", c) + if suffix: + c = c[: -len(suffix[0])] + return _validate_curie(c, prefixes) + + +def _validate_column(column, column_name, df_name, schema_def): + """Given a schema definition and the column of a dataframe, verify that the column satifies + the schema. + """ + + errors = [] + + if schema_def.get("unique"): + if column.nunique() != len(column): + errors.append(f"Column {column_name} in dataframe {df_name} is not unique.") + + if "nullable" in schema_def and not schema_def["nullable"]: + if any(_is_null(v) for v in column): + errors.append(f"Column {column_name} in dataframe {df_name} contains empty values.") + + if schema_def.get("type") == "human-readable string": + non_readables = [v for v in column if not _validate_human_readable_string(v)] + if non_readables: + errors.append( + f"Column {column_name} in dataframe {df_name} contains non-human-readable " + f"values like {non_readables[0]}" + ) + + if schema_def.get("type") in ("curie", "suffixed curie"): + validation_func = _validate_curie if schema_def.get("type") == "curie" else _validate_suffixed_curie + non_valid_curies = [v for v in column if not validation_func(v, schema_def.get("prefixes"))] + if non_valid_curies: + errors.append( + f"Column {column_name} in dataframe {df_name} contains invalid ontology values like " + f"{non_valid_curies[0]}." + ) + if "prefixes" in schema_def: + errors[-1] += f" Values must be curies from one of these ontologies {schema_def['prefixes']}." + + if "enum" in schema_def: + bad_enums = [v for v in column if v not in schema_def["enum"]] + if bad_enums: + errors.append( + f"Column {column_name} in dataframe {df_name} contains unpermitted values like " + f"{bad_enums[0]}. Values must be one of {schema_def['enum']}." + ) + + return errors + + +def _validate_dict(dict_, dict_name, schema_def): + """Given a schema definition and dict, verify that the dict satifies the schema.""" + + errors = [] + + for key in schema_def.get("keys", []): + if key not in dict_: + errors.append(f"{dict_name} is missing key {key}.") + elif schema_def["keys"][key]: + if schema_def["keys"][key]["type"] == "stringified list of dicts": + if not _validate_stringified_list_of_dicts(dict_[key]): + errors.append( + f"Key {key} in {dict_name} should be a JSON-encoded list of dicts, but it is {dict_[key]}" + ) + elif schema_def["keys"][key]["type"] == "dict": + errors.extend(_validate_dict(dict_[key], key, schema_def["keys"][key])) + elif schema_def["keys"][key]["type"] == "curie": + if not _validate_curie(dict_[key], schema_def["keys"][key]["prefixes"]): + errors.append(f"Key {key} in {dict_name} contains invalid ontology value.") + if "nullable" in schema_def["keys"][key] and not schema_def["keys"][key]["nullable"]: + if _is_null(dict_[key]): + errors.append(f"Key {key} in dict {dict_name} is an empty value.") + + return errors + + +def _validate_dataframe(df, df_name, schema_def): + """Given a dataframe and schema definition, verify that the dataframe follows the schema.""" + + errors = [] + + if "index" in schema_def: + errors.extend(_validate_column(df.index, "index", df_name, schema_def["index"])) + + for column in schema_def.get("columns", []): + if column not in df.columns: + errors.append(f"Dataframe {df_name} is missing column {column}.") + else: + errors.extend(_validate_column(df[column], column, df_name, schema_def["columns"][column])) + + return errors + + +def get_schema_definition(version): + """Look up and read a schema definition based on a version number like "1.0.0".""" + + path = os.path.join( + os.path.dirname(os.path.realpath(__file__)), "schema_definitions", version.replace(".", "_") + ".yaml" + ) + + if not os.path.isfile(path): + raise ValueError(f"No definition for version {version} found.") + + return yaml.load(open(path), Loader=yaml.FullLoader) + + +def deep_check(adata, schema_def): + """Perform a "deep" check of the AnnData object using the schema definition. + + This checks all the columns and unstructured metadata rather than just the version. + + Returns a list of error messages. If that list is empty, the object passed validation. + """ + + errors = [] + + for component, component_def in schema_def["components"].items(): + if component_def["type"] == "dataframe": + errors.extend(_validate_dataframe(getattr(adata, component), component, component_def)) + elif component_def["type"] == "dict": + errors.extend(_validate_dict(getattr(adata, component), component, component_def)) + else: + raise ValueError(f"Unexpected component type {component['type']}") + + return errors + + +def validate_adata(adata, shallow): + """Validate an AnnData object. If shallow, just check that the required version information is + present. + """ + + # Does it have the version information written into uns? + if "version" not in adata.uns_keys() or "corpora_schema_version" not in adata.uns["version"]: + print("AnnData file is missing corpora version information") + return False + + # We can stop here if it's a "shallow" check, that is, if we're just + # checking that version is present. + if shallow: + return True + + schema_def = get_schema_definition(adata.uns["version"]["corpora_schema_version"]) + + errors = deep_check(adata, schema_def) + + for error in errors: + print(error) + + return not errors + + +def validate(h5ad_path, shallow=False): + """Entry point for validation.""" + + try: + import scanpy + except ImportError: + raise ImportError("scanpy must be installed for cellxgene schema") + + try: + adata = scanpy.read_h5ad(h5ad_path, backed="r") + except (OSError, TypeError): + print(f"Unable to open {h5ad_path} with scanpy.") + sys.exit(1) + + if not validate_adata(adata, shallow): + sys.exit(1) diff --git a/server/test/fixtures/hgnc_example.txt.gz b/server/test/fixtures/hgnc_example.txt.gz new file mode 100644 index 00000000..4c7704c4 Binary files /dev/null and b/server/test/fixtures/hgnc_example.txt.gz differ diff --git a/server/test/fixtures/schema_test_data/generate_test_data.sh b/server/test/fixtures/schema_test_data/generate_test_data.sh new file mode 100755 index 00000000..d0f7d699 --- /dev/null +++ b/server/test/fixtures/schema_test_data/generate_test_data.sh @@ -0,0 +1,139 @@ +#!/bin/bash +wget "https://s3-us-west-2.amazonaws.com/10x.files/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz" +tar xf "pbmc3k_filtered_gene_bc_matrices.tar.gz" + +python3 - < genes_tmp.tsv; mv genes_tmp.tsv merged/genes.tsv + +echo -e "\n\n\nRunning tutorial on original\n\n\n" +Rscript - <