diff --git a/requirements.txt b/requirements.txt index 42390d9e9..4399e0823 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,10 +2,9 @@ cloudpathlib[all]>=0.16.0 cyvcf2>=0.30.18 dill>=0.3.7 hail==0.2.133 +hpo3 httpx>=0.27.0 Jinja2>=3.1.3 -networkx>=3 -obonet>=1 pandas>=2 peds>=1.2.0 phenopackets>=2 @@ -13,7 +12,5 @@ protobuf==3.20.2 pydantic>=2.5.2 pyspark>=3.5.1 python-dateutil>=2 -# this is a ruff package with a python interface, limited to specific python versions -semsimian tabulate>=0.8.9 toml==0.10.2 diff --git a/src/talos/CPG/MakePhenopackets.py b/src/talos/CPG/MakePhenopackets.py index b54dc5362..03259a2ae 100755 --- a/src/talos/CPG/MakePhenopackets.py +++ b/src/talos/CPG/MakePhenopackets.py @@ -15,17 +15,17 @@ import re from argparse import ArgumentParser +from collections import defaultdict -import networkx as nx import phenopackets.schema.v2 as pps2 from google.protobuf.json_format import MessageToJson -from obonet import read_obo +from pyhpo import Ontology from metamist.graphql import gql, query from talos.static_values import get_logger HPO_KEY = 'HPO Terms (present)' -HPO_RE = re.compile(r'HP:[0-9]+') +HPO_RE = re.compile(r'HP:\d+') PARTICIPANT_QUERY = gql( """ query MyQuery($project: String!, $sequencing_type: String!, $technology: String!) { @@ -48,11 +48,14 @@ }""", ) +# create an ontology object, once +_ = Ontology() + # map the integer reported sex values to the enum reported_sex_map = {1: pps2.Sex.MALE, 2: pps2.Sex.FEMALE} -def find_hpo_labels(metamist_data: dict, hpo_file: str | None = None) -> dict[str, list[dict[str, str]]]: +def find_hpo_labels(metamist_data: dict) -> dict[str, list[dict[str, str]]]: """ match HPO terms to their plaintext names @@ -62,46 +65,34 @@ def find_hpo_labels(metamist_data: dict, hpo_file: str | None = None) -> dict[st disease genes and their MOI terms directly, which messes with what we are trying to do: associate patient _phenotypes_ with variant _genes_. + As we're now using hpo3/pyHPO which ships with a built-in ontology, we can use assume presence of a valid Ontology, + simplifying the code + Args: metamist_data (): - hpo_file (): Returns: dict, participant IDs to HPO:labels """ - all_hpos: set[str] = set() - per_sg_hpos: dict[str, set[str]] = {} - - moi_nodes: set[str] = set() - hpo_graph = None - if hpo_file: - # create a graph of HPO terms - hpo_graph = read_obo(hpo_file, ignore_obsolete=False) - moi_nodes = nx.ancestors(hpo_graph, 'HP:0000005') + per_sg_hpos: dict[str, list[dict[str, str]]] = defaultdict(list) for sg in metamist_data['project']['sequencingGroups']: - hpos = set(HPO_RE.findall(sg['sample']['participant']['phenotypes'].get(HPO_KEY, ''))) - - # groom out any strictly MOI related terms - hpos -= moi_nodes - - all_hpos.update(hpos) - per_sg_hpos[sg['id']] = hpos - - # no label obtained, that's fine... - if not (hpo_file and hpo_graph): - return {sg: [{'id': hp, 'label': 'Unknown'} for hp in hpos] for sg, hpos in per_sg_hpos.items()} + # select all HPO terms, so long as they are not a child of 'Mode of Inheritance' (HP:0000005) + hpos = { + hpo_term + for hpo_term in HPO_RE.findall(sg['sample']['participant']['phenotypes'].get(HPO_KEY, '')) + if 'HP:0000005' not in Ontology.get_hpo_object(hpo_term).all_parents + } - # create a dictionary of HPO terms to their text - hpo_to_text: dict[str, str] = {} - for hpo in all_hpos: - if not hpo_graph.has_node(hpo): - get_logger(__file__).error(f'HPO term was absent from the tree: {hpo}') - hpo_to_text[hpo] = 'Unknown' - else: - hpo_to_text[hpo] = hpo_graph.nodes[hpo]['name'] + # allow for HPO terms to be missing from this edition of the ontology + for hpo in hpos: + try: + per_sg_hpos[sg['id']].append({'id': hpo, 'label': Ontology.get_hpo_object(hpo).name}) + except ValueError: + get_logger(__file__).error(f'HPO term was absent from the tree: {hpo}') + per_sg_hpos[sg['id']].append({'id': hpo, 'label': 'Unknown'}) - return {sg: [{'id': hp, 'label': hpo_to_text[hp]} for hp in hpos] for sg, hpos in per_sg_hpos.items()} + return per_sg_hpos def assemble_phenopackets(dataset: str, metamist_data: dict, hpo_lookup: dict[str, list[dict[str, str]]]): @@ -170,7 +161,7 @@ def assemble_phenopackets(dataset: str, metamist_data: dict, hpo_lookup: dict[st return cohort -def main(output: str, dataset: str, seq_type: str, tech: str = 'short-read', hpo_file: str | None = None): +def main(output: str, dataset: str, seq_type: str, tech: str = 'short-read'): """ Assemble a cohort phenopacket from the metamist data Args: @@ -178,7 +169,6 @@ def main(output: str, dataset: str, seq_type: str, tech: str = 'short-read', hpo dataset (str): the dataset to query for seq_type (str): exome/genome tech (str): type of sequence data to query for - hpo_file (str): path to an obo file containing HPO tree """ # pull all the relevant data from metamist @@ -188,7 +178,7 @@ def main(output: str, dataset: str, seq_type: str, tech: str = 'short-read', hpo ) # match names to HPO terms - labelled_hpos = find_hpo_labels(hpo_file=hpo_file, metamist_data=metamist_data) + labelled_hpos = find_hpo_labels(metamist_data=metamist_data) # build the cohort cohort = assemble_phenopackets(dataset=dataset, metamist_data=metamist_data, hpo_lookup=labelled_hpos) @@ -203,11 +193,10 @@ def cli_main(): parser.add_argument('--dataset', help='The dataset to query for') parser.add_argument('--output', help='The output file') parser.add_argument('--type', help='Sequencing type (exome or genome)') - parser.add_argument('--hpo', help='HPO ontology file') parser.add_argument('--tech', help='Sequencing technology', default='short-read') args = parser.parse_args() - main(dataset=args.dataset, output=args.output, seq_type=args.type, hpo_file=args.hpo, tech=args.tech) + main(dataset=args.dataset, output=args.output, seq_type=args.type, tech=args.tech) if __name__ == '__main__': diff --git a/src/talos/FindGeneSymbolMap.py b/src/talos/FindGeneSymbolMap.py index 01d4d1503..be0a6df1b 100644 --- a/src/talos/FindGeneSymbolMap.py +++ b/src/talos/FindGeneSymbolMap.py @@ -34,7 +34,8 @@ async def match_ensgs_to_symbols(genes: list[str], session: ClientSession) -> di r.raise_for_status() json_reponse = await r.json() # match symbol to the ENSG (or Unknown if the key is missing, or has a None value) - return {value.get('display_name'): key for key, value in json_reponse.items() if value} + # we want the key to be ENSG + return {key: value.get('display_name') for key, value in json_reponse.items() if value} async def match_symbol_to_ensg(gene_symbol: str, session: ClientSession) -> tuple[str, str]: diff --git a/src/talos/GeneratePanelData.py b/src/talos/GeneratePanelData.py index 7807e0f7a..4da414816 100644 --- a/src/talos/GeneratePanelData.py +++ b/src/talos/GeneratePanelData.py @@ -12,8 +12,7 @@ import phenopackets.schema.v2 as pps2 from google.protobuf.json_format import ParseDict -from networkx import MultiDiGraph -from obonet import read_obo +from pyhpo import Ontology from talos.config import config_retrieve from talos.models import ParticipantHPOPanels, PhenoPacketHpo, PhenotypeMatchedPanels @@ -33,6 +32,9 @@ PANELS_ENDPOINT = PANELAPP_HARD_CODED_DEFAULT DEFAULT_PANEL = PANELAPP_HARD_CODED_BASE_PANEL +# create an ontology object, once +_ = Ontology() + def get_panels(endpoint: str = PANELS_ENDPOINT) -> dict[str, set[int]]: """ @@ -89,60 +91,15 @@ def set_up_cohort_pmp(cohort: pps2.Cohort) -> tuple[PhenotypeMatchedPanels, set[ return hpo_dict, all_hpos -def match_hpo_terms( - panel_map: dict[str, set[int]], - hpo_tree: MultiDiGraph, - hpo_str: str, - selections: set[int] | None = None, -) -> set[int]: - """ - get panels relevant for this HPO using a recursive edge traversal - for live terms we recurse on all parents - if a term is obsolete we instead check each replacement term - - relevant usage guide: - https://github.com/dhimmel/obonet/blob/main/examples/go-obonet.ipynb - - Args: - panel_map (dict): - hpo_tree (): a graph object representing the HPO tree - hpo_str (str): the query HPO term - selections (set[int]): collected panel IDs so far - - Returns: - set: panel IDs relating to this HPO term - """ - - if selections is None: - selections = set() - - # identify identical match and select the panel - if hpo_str in panel_map: - selections.update(panel_map[hpo_str]) - - # if a node is invalid, recursively call this method for each replacement D: - # there are simpler ways, just none that are as fun to write - if not hpo_tree.has_node(hpo_str): - get_logger().error(f'HPO term was absent from the tree: {hpo_str}') - return selections - - hpo_node = hpo_tree.nodes[hpo_str] - if hpo_node.get('is_obsolete', 'false') == 'true': - for hpo_term in hpo_node.get('replaced_by', []): - selections.update(match_hpo_terms(panel_map, hpo_tree, hpo_term, selections)) - # search for parent(s), even if the term is obsolete - for hpo_term in hpo_node.get('is_a', []): - selections.update(match_hpo_terms(panel_map, hpo_tree, hpo_term, selections)) - return selections - - -def match_hpos_to_panels(hpo_panel_map: dict[str, set[int]], hpo_file: str, all_hpos: set[str]) -> dict[str, set[int]]: +def match_hpos_to_panels(hpo_panel_map: dict[str, set[int]], all_hpos: set[str]) -> dict[str, set[int]]: """ take the HPO terms from the participant metadata, and match to panels + greatly simplified implementation, as we're now using hpo3/pyHPO which ships with a built-in ontology + we're checking for clashes between PanelApp contents and each exact term, or its ancestors + Args: hpo_panel_map (dict): panel IDs to all related panels - hpo_file (str): path to an obo file containing HPO tree all_hpos (set[str]): collection of all unique HPO terms Returns: @@ -150,12 +107,23 @@ def match_hpos_to_panels(hpo_panel_map: dict[str, set[int]], hpo_file: str, all_ a second dictionary linking all HPO terms to their plaintext names """ - hpo_graph = read_obo(hpo_file, ignore_obsolete=False) - - hpo_to_panels = {} + hpo_to_panels: dict[str, set[int]] = {} for hpo in all_hpos: - panel_ids = match_hpo_terms(panel_map=hpo_panel_map, hpo_tree=hpo_graph, hpo_str=hpo) - hpo_to_panels[hpo] = panel_ids + # set must exist, even if empyt - defaultdict doesn't do this + hpo_to_panels[hpo] = set() + + # if we have an exact match, use that as a starting point + if hpo in hpo_panel_map: + hpo_to_panels[hpo] = hpo_panel_map[hpo] + + # try and fetch the node in the ontology + try: + hpo_node = Ontology.get_hpo_object(hpo) + for parent in hpo_node.all_parents: + if parent.id in hpo_panel_map: + hpo_to_panels[hpo].update(hpo_panel_map[parent.id]) + except ValueError: + get_logger(__file__).error(f'HPO term was absent from the tree: {hpo}') return hpo_to_panels @@ -187,12 +155,11 @@ def cli_main(): parser = ArgumentParser() parser.add_argument('--input', help='GA4GH Cohort/PhenoPacket File') parser.add_argument('--output', help='Path to write PhenotypeMatchedPanels to (JSON)') - parser.add_argument('--hpo', help='Local copy of HPO obo file', required=True) args = parser.parse_args() - main(ga4gh_cohort_file=args.input, panel_out=args.output, hpo_file=args.hpo) + main(ga4gh_cohort_file=args.input, panel_out=args.output) -def main(ga4gh_cohort_file: str, panel_out: str, hpo_file: str): +def main(ga4gh_cohort_file: str, panel_out: str): """ query PanelApp - get all panels and their assc. HPOs read Cohort/PhenoPacket file @@ -202,7 +169,6 @@ def main(ga4gh_cohort_file: str, panel_out: str, hpo_file: str): Args: ga4gh_cohort_file (str): path to GA4GH Cohort/PhenoPacket file panel_out (str): where to write PhenotypeMatchedPanels file - hpo_file (str): path to an obo file containing HPO tree """ # read the Cohort/PhenoPacket file as a JSON @@ -217,7 +183,7 @@ def main(ga4gh_cohort_file: str, panel_out: str, hpo_file: str): # match HPO terms to panel IDs # returns a lookup of each HPO term in the cohort, and panels it is associated with - hpo_to_panels = match_hpos_to_panels(hpo_panel_map=panels_by_hpo, all_hpos=all_hpos, hpo_file=hpo_file) + hpo_to_panels = match_hpos_to_panels(hpo_panel_map=panels_by_hpo, all_hpos=all_hpos) match_participants_to_panels(pmp_dict, hpo_to_panels) diff --git a/src/talos/HPOFlagging.py b/src/talos/HPOFlagging.py index 4acb904f5..d29770c46 100644 --- a/src/talos/HPOFlagging.py +++ b/src/talos/HPOFlagging.py @@ -1,14 +1,18 @@ """ Takes the result date from ValidateMOI -For each reportable result, identify if the variant is a strong phenotype match for the family +For each reportable result, identify if the variant gene is enriched relateive to the participant's phenotypes +This uses the HPO3/pyHPO gene enrichment model, see: -1. Write 2 results - the result set annotated with Phenotype-match flags -2. The annotated result set, reduced to only phenotype matches +https://github.com/anergictcell/pyhpo?tab=readme-ov-file#get-genes-enriched-in-an-hposet +https://pyhpo.readthedocs.io/en/latest/stats.html#statistics + +This pretty cleanly matches our use-case: "You have a set of HPOTerms and want to find the most likely causative gene" -Include something here about where to find the inputs ... +Here we're doing an enrichment test using the participant's phenotypes, then accepting any genes with an enrichment +probability below 0.05 -gene_to_phtnotype: Download: https://hpo.jax.org/app/data/annotations#:~:text=download-,GENES,-TO%20PHENOTYPE" -phenio.db: Download: https://data.monarchinitiative.org/monarch-kg/latest/phenio.db.gz +1. Write 2 results - the result set annotated with Phenotype-match flags +2. The annotated result set, reduced to only phenotype matches The removal of variants where a phenotype match is required but not found it now done here A list of categories in config is used to determine which variants are required to have a phenotype match @@ -19,139 +23,89 @@ """ from argparse import ArgumentParser -from collections import defaultdict -from semsimian import Semsimian +from pyhpo import Gene, HPOSet, Ontology +from pyhpo.stats import EnrichmentModel from talos.config import config_retrieve from talos.models import ParticipantResults, ResultData from talos.static_values import get_granular_date -from talos.utils import phenotype_label_history, read_json_from_path +from talos.utils import get_logger, phenotype_label_history, read_json_from_path -_SEMSIM_CLIENT: Semsimian | None = None +# some setup for hpo3 +_ = Ontology() +ENRICHMENT_MODEL = EnrichmentModel('gene') -def get_sem_client(phenio_db: str | None = None) -> Semsimian: - """ - create or retrieve a Semsimian client - - Args: - phenio_db (str | None): needs to be present for the first call - Returns: - Semsimian instance +def get_gene_enrichment(hpo_terms: list[str]) -> dict[str, dict[str, float]]: """ - global _SEMSIM_CLIENT - if _SEMSIM_CLIENT is None: - _SEMSIM_CLIENT = Semsimian(spo=None, predicates=['rdfs:subClassOf'], resource_path=phenio_db) - return _SEMSIM_CLIENT - - -def parse_genes_to_hpo(g2p_file: str, ensg2symbol: str, genes: set[str]) -> dict[str, set[str]]: - """ - Parse genes to phenotype file from Jax. - Returns a dict of gene_symbol -> set of HPO ids - + get the gene enrichment for a set of HPO terms Args: - g2p_file (): - ensg2symbol (): the gene symbols to Ensembl gene IDs for genes relevant to this analysis - genes (set[str]): all genes in this report + hpo_terms (): Returns: - dict[str, set[str]]: ENSG -> set of HPO ids + dictionary of genes and their enrichments """ - # lookup of {symbol: ENSG} - ensg_map: dict[str, str] = read_json_from_path(ensg2symbol) + # pull the enrichment threshold from config + enrichment_threshold = config_retrieve(['HPOFlagging', 'enrichment_threshold'], 0.05) - gene_to_phenotype: dict[str, set[str]] = defaultdict(set) - with open(g2p_file, encoding='utf-8') as f: - for line in f: - ncbi_gene_id, gene_symbol, hpo_id, hpo_name, frequency, disease_id = line.split('\t') - if gene_symbol in ensg_map and ensg_map[gene_symbol] in genes: - gene_to_phenotype.setdefault(ensg_map[gene_symbol], set()).add(hpo_id) - return gene_to_phenotype + # check for genes enriched in the participant's HPO terms + participant_hpo3 = HPOSet.from_queries(hpo_terms) + # find genes which are enriched relative to the participant's HPO terms + filtered_enriched: dict[str, dict[str, float]] = { + gene['item'].name: {'enrichment': gene['enrichment'], 'fold': gene['fold']} + for gene in ENRICHMENT_MODEL.enrichment(method='hypergeom', hposet=participant_hpo3) + if gene['enrichment'] <= enrichment_threshold + } -def find_genes_in_these_results(result_object: ResultData) -> set[str]: - """ - identify all the genes we care about finding phenotype matches for + return filtered_enriched - Args: - result_object (ResultData): - Returns: - set[str], all ensembl gene IDs in - """ - ensgs: set[str] = set() - for participant in result_object.results.values(): - for variant in participant.variants: - # get the gene ID - ensgs.add(variant.var_data.info['gene_id']) # type: ignore[arg-type] - - # for structural variants, add all the LOF'd genes - # TODO (mwelland): if/when we create other SV categories, we may need to catch those here too - if lof := variant.var_data.info.get('predicted_lof'): - ensgs.update(set(str(lof).split(','))) - - return ensgs - - -def annotate_phenotype_matches(result_object: ResultData, gen_phen: dict[str, set[str]]): +def annotate_phenotype_matches(result_object: ResultData, ensg2symbol: str): """ - for each variant, find any phenotype matches between the participant and gene HPO sets - - checks config for HPOFlagging.strict (bool) - if strict, this will test for an exact overlap - if strict is False, this will be done semantically via Semsimian + for each variant, find probability of the variant gene showing an enriched association with the + set of patient phenotype terms Args: result_object (ResultData): - gen_phen (dict): mapping of ENSGs to relevant HPO terms + ensg2symbol (str): path to the ENSG to symbol mapping """ - semantic_match = config_retrieve(['HPOFlagging', 'semantic_match'], False) - - min_similarity: float = config_retrieve(['HPOFlagging', 'min_similarity']) + ensg_map: dict[str, str] = read_json_from_path(ensg2symbol) for participant in result_object.results.values(): + # no HPOs, no problems + if not participant.metadata.phenotypes: + continue + participant_hpos_dict = {hpo.id: hpo.label for hpo in participant.metadata.phenotypes} - participant_hpos = set(participant_hpos_dict.keys()) + participant_hpos = list(participant_hpos_dict.keys()) + + # find genes which are enriched relative to the participant's HPO terms + filtered_enriched = get_gene_enrichment(participant_hpos) for variant in participant.variants: - var_gene = variant.var_data.info['gene_id'] - gene_hpos = gen_phen.get(var_gene, set()) # type: ignore[arg-type] - - # under strict matching we require exact overlapping terms - # we always run a strict match - for hpo_id in participant_hpos & gene_hpos: - variant.phenotype_labels.add(f'{hpo_id}: {participant_hpos_dict[hpo_id]}') - - # optionally also use semantic matching for phenotypic similarity - if participant_hpos and gene_hpos and semantic_match: - termset_similarity = get_sem_client().termset_pairwise_similarity(participant_hpos, gene_hpos) - # Convert object terms (gene_phenotypes) to lookup dict - object_termset = { - term_dict['id']: term_dict['label'] - for term in termset_similarity['object_termset'] - for term_dict in term.values() - } - - # Find all phenotype matches that meet the min_score threshold - pheno_matches = { - f"{match['object_id']}: {object_termset[match['object_id']]}" - for match in termset_similarity['subject_best_matches']['similarity'].values() - if float(match['ancestor_information_content']) > min_similarity - } - - # skip if no matches - don't assign a date if there are no matches - if not pheno_matches: - continue + var_symbol = ensg_map.get(variant.var_data.info['gene_id']) # type: ignore[assignment] + if var_symbol in filtered_enriched: + variant.phenotype_labels.add(f'p.{filtered_enriched[var_symbol]["enrichment"]:.4f}') if variant.date_of_phenotype_match is None: variant.date_of_phenotype_match = get_granular_date() - variant.phenotype_labels = pheno_matches + # get exactly matching HPO terms for the gene + try: + gene_hpo_set = Gene.get(var_symbol).hpo_set() + gene_hpo_strings = {hpo_term.id for hpo_term in gene_hpo_set} + for hpo_term in gene_hpo_strings & set(participant_hpos): + if variant.date_of_phenotype_match is None: + variant.date_of_phenotype_match = get_granular_date() + hpo_object = Ontology.get_hpo_object(hpo_term) + variant.phenotype_labels.add(f'{hpo_term}: {hpo_object.name} (exact)') + except KeyError: + get_logger(__file__).warning(f'Could not find gene {var_symbol} in the ontology') phenotype_label_history(result_object) @@ -225,16 +179,12 @@ def cli_main(): parser = ArgumentParser(description='') parser.add_argument('--input', help='The Result data in JSON form') parser.add_argument('--gene_map', help='A map of gene symbol to ENSG for all genes in this analysis') - parser.add_argument('--gen2phen', help='path to the genotype-phenotype file') - parser.add_argument('--phenio', help='A phenio DB file') parser.add_argument('--output', help='Annotated full output') parser.add_argument('--phenout', help='Annotated phenotype-only output') args = parser.parse_args() main( result_file=args.input, gene_map=args.gene_map, - gen2phen=args.gen2phen, - phenio=args.phenio, out_path=args.output, phenout=args.phenout, ) @@ -243,8 +193,6 @@ def cli_main(): def main( result_file: str, gene_map: str, - gen2phen: str, - phenio: str, out_path: str, phenout: str | None = None, ): @@ -253,8 +201,6 @@ def main( Args: result_file (str): path to the ValidateMOI output JSON gene_map (str): output of FindGeneSymbolMap, JSON - gen2phen (str): path to a test file of known Phenotypes per gene - phenio (str): path to a PhenoIO DB file out_path (str): path to write the annotated results phenout (str): optional, path to phenotype filtered outputs """ @@ -262,19 +208,8 @@ def main( # read the results JSON into an object results = read_json_from_path(result_file, return_model=ResultData) - # find all the genes present in this report - relevant_ensgs = find_genes_in_these_results(results) - - # for each gene, identify the HPO terms that are associated with it - gen_phen_dict = parse_genes_to_hpo(gen2phen, gene_map, genes=relevant_ensgs) - - # unless strict, set up the client with a phenio file - # this is a bit hacky, but means we don't need to pass the phenio path around - if not config_retrieve(['HPOFlagging', 'strict'], False): - _semsim = get_sem_client(phenio_db=phenio) - # label phenotype matches, and write the results to a JSON file - annotate_phenotype_matches(results, gen_phen_dict) + annotate_phenotype_matches(results, gene_map) # remove any variants where a phenotype match is required but not found remove_phenotype_required_variants(results) diff --git a/src/talos/models.py b/src/talos/models.py index 2fe94b621..93aafe3a5 100644 --- a/src/talos/models.py +++ b/src/talos/models.py @@ -95,7 +95,7 @@ class VariantCommon(BaseModel): boolean_categories: list[str] = Field(default_factory=list, exclude=True) sample_categories: list[str] = Field(default_factory=list, exclude=True) sample_support: list[str] = Field(default_factory=list, exclude=True) - phased: dict = Field(default_factory=dict) + phased: dict = Field(default_factory=dict, exclude=True) def __str__(self): return repr(self) diff --git a/test/conftest.py b/test/conftest.py index 898be00bb..2db1aa7d1 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -10,6 +10,7 @@ import hail as hl import pytest from cyvcf2 import VCFReader +from pyhpo import Ontology from talos.data_model import BaseFields, Entry, SneakyTable, TXFields, VepVariant from talos.utils import create_small_variant, read_json_from_path @@ -24,7 +25,6 @@ LABELLED = join(INPUT, '1_labelled_variant.vcf.bgz') Talos_OUTPUT = join(INPUT, 'aip_output_example.json') DE_NOVO_PED = join(INPUT, 'de_novo_ped.fam') -FAKE_OBO = join(INPUT, 'hpo_test.obo') LOOKUP_PED = join(INPUT, 'mock_sm_lookup.json') PHASED_TRIO = join(INPUT, 'newphase.vcf.bgz') PED_FILE = join(INPUT, 'pedfile.ped') @@ -59,6 +59,17 @@ def fixture_hail_cleanup(): filename.unlink() +@pytest.fixture(name='setup_ontology', scope='session', autouse=True) +def fixture_set_up_pyhpo_ontology(): + """ + a fixture to instantiate a pyhpo ontology, once + """ + _ = Ontology() + + # yield something to suspend + yield '' + + @pytest.fixture(name='make_a_mt', scope='session') def fixture_make_a_mt(tmp_path_factory) -> hl.MatrixTable: """ @@ -99,12 +110,6 @@ def fixture_test_input_models_path() -> str: return join(INPUT, 'models') -@pytest.fixture(name='fake_obo_path', scope='session') -def fixture_fake_obo() -> str: - """path to fake obo""" - return FAKE_OBO - - @pytest.fixture(name='fake_panelapp_overview', scope='session') def fixture_panelapp_overview() -> Any: """path to panelapp_overview json""" diff --git a/test/input/hpo_test.obo b/test/input/hpo_test.obo deleted file mode 100644 index 7466a28e9..000000000 --- a/test/input/hpo_test.obo +++ /dev/null @@ -1,49 +0,0 @@ -format-version: 1.2 -ontology: test.ob - -[Term] -id: HP:1 -name: Philosopher's Stone -comment: root node - -[Term] -id: HP:2 -name: Chamber of Secrets -comment: first sub_node -is_a: HP:1 - -[Term] -id: HP:3 -name: Prisoner of Azkaban -comment: where my Hippogriff at? -is_a: HP:2 - -[Term] -id: HP:4 -name: Goblet of Fire -comment: https://www.youtube.com/watch?v=2g3xrYjreN4 -is_a: HP:3 - -[Term] -id: HP:5 -name: Order of the phoenix -comment: https://www.youtube.com/watch?v=5JqY-6q-RNA -is_a: HP:4 - -[Term] -id: HP:6 -name: Half-blood Prince -comment: TL;DR it was Snape -is_a: HP:5 - -[Term] -id: HP:7a -name: Deathly Hallows -comment: the one with the cool cartoon intro -is_a: HP:6 - -[Term] -id: HP:7b -name: Deathly Hallows 2 -comment: it's all ogre now -is_a: HP:6 diff --git a/test/test_hpo_flagging.py b/test/test_hpo_flagging.py deleted file mode 100644 index 45a359e6d..000000000 --- a/test/test_hpo_flagging.py +++ /dev/null @@ -1,26 +0,0 @@ -from talos.HPOFlagging import find_genes_in_these_results -from talos.models import Coordinates, ParticipantMeta, ParticipantResults, ReportVariant, ResultData, SmallVariant - -TEST_COORDS = Coordinates(chrom='1', pos=1, ref='A', alt='C') -SMALL_1 = SmallVariant( - info={'gene_id': 'ENSG1'}, - coordinates=TEST_COORDS, - transcript_consequences=[], -) -SV_1 = SmallVariant( - info={'gene_id': 'ENSG2', 'predicted_lof': 'ENSG3,ENSG4'}, - coordinates=TEST_COORDS, - transcript_consequences=[], -) - - -def test_find_genes_in_these_results(): - pm = ParticipantMeta(ext_id='male', family_id='family_1') - test_results = ResultData( - results={ - 'male': ParticipantResults(metadata=pm, variants=[ReportVariant(sample='male', var_data=SMALL_1)]), - 'female': ParticipantResults(metadata=pm, variants=[ReportVariant(sample='male', var_data=SV_1)]), - }, - ) - identified = find_genes_in_these_results(test_results) - assert identified == {'ENSG1', 'ENSG2', 'ENSG3', 'ENSG4'} diff --git a/test/test_metamist_hpo.py b/test/test_metamist_hpo.py index 99725f465..7564ccf95 100644 --- a/test/test_metamist_hpo.py +++ b/test/test_metamist_hpo.py @@ -2,10 +2,7 @@ test file for metamist panel-participant matching """ -import networkx as nx -from obonet import read_obo - -from talos.GeneratePanelData import get_panels, match_hpo_terms, match_hpos_to_panels, match_participants_to_panels +from talos.GeneratePanelData import get_panels, match_hpos_to_panels, match_participants_to_panels from talos.models import ParticipantHPOPanels, PhenotypeMatchedPanels @@ -18,54 +15,14 @@ def test_get_panels(httpx_mock, fake_panelapp_overview): assert panels_parsed == {'HP:1': {2}, 'HP:4': {1}, 'HP:6': {2}} -def test_match_hpo_terms(fake_obo_path): - """ - check that HP tree traversal works - this test is kinda limited now that the layer count is constant - """ - obo_parsed = read_obo(fake_obo_path) - panel_map = {'HP:2': {1, 2}} - assert match_hpo_terms(panel_map=panel_map, hpo_tree=obo_parsed, hpo_str='HP:4') == {1, 2} - assert match_hpo_terms(panel_map=panel_map, hpo_tree=obo_parsed, hpo_str='HP:2') == {1, 2} - - -def test_match_hpos_to_panels(fake_obo_path): +def test_match_hpos_to_panels(): """ test the hpo-to-panel matching - this has now been adjusted to account for the full leaf-to-root traversal, instead of stopping at 3 layers - """ - panel_map = {'HP:2': {1, 2}, 'HP:5': {5}} - assert match_hpos_to_panels(panel_map, fake_obo_path, all_hpos={'HP:4', 'HP:7a'}) == { - 'HP:4': {1, 2}, - 'HP:7a': {1, 2, 5}, - } - - # full depth from the terminal node should capture all panels - assert match_hpos_to_panels(panel_map, fake_obo_path, all_hpos={'HP:4', 'HP:7a'}) == { - 'HP:4': {1, 2}, - 'HP:7a': {1, 2, 5}, - } - - -def test_read_hpo_tree(fake_obo_path): - """ - check that reading the obo tree works """ - obo_parsed = read_obo(fake_obo_path) - assert isinstance(obo_parsed, nx.MultiDiGraph) - assert list(nx.bfs_edges(obo_parsed, 'HP:1', reverse=True)) == [ - ('HP:1', 'HP:2'), - ('HP:2', 'HP:3'), - ('HP:3', 'HP:4'), - ('HP:4', 'HP:5'), - ('HP:5', 'HP:6'), - ('HP:6', 'HP:7a'), - ('HP:6', 'HP:7b'), - ] - assert obo_parsed.nodes()['HP:3'] == { - 'name': 'Prisoner of Azkaban', - 'comment': 'where my Hippogriff at?', - 'is_a': ['HP:2'], + panel_map = {'HP:0000574': {1, 2}, 'HP:0000234': {5}} + assert match_hpos_to_panels(panel_map, all_hpos={'HP:0000574', 'HP:0000234'}) == { + 'HP:0000574': {1, 2, 5}, + 'HP:0000234': {5}, }