diff --git a/clickhouse_search/managers.py b/clickhouse_search/managers.py index af5d5c4d28..b514d2620d 100644 --- a/clickhouse_search/managers.py +++ b/clickhouse_search/managers.py @@ -22,6 +22,10 @@ from seqr.utils.xpos_utils import get_xpos, MIN_POS, MAX_POS, CHROMOSOME_CHOICES +class InvalidDatasetTypeException(Exception): + pass + + class SearchQuerySet(QuerySet): @property @@ -414,37 +418,30 @@ def populations(self): def _filter_frequency(self, results, **kwargs): return results - def _filter_annotations(self, results, annotations=None, pathogenicity=None, genes=None, intervals=None, exclude_locations=False, **kwargs): + def _filter_annotations(self, results, annotations_filter_q=None, transcript_filters=None, genes=None, intervals=None, exclude_locations=False, **kwargs): if exclude_locations and genes: intervals = self._format_gene_intervals(genes) genes = None - results = self._filter_locations(results, genes, intervals, exclude_locations=exclude_locations, **kwargs) - filter_qs, transcript_filters = self._parse_annotation_filters(annotations or {}, pathogenicity) if (annotations or pathogenicity) else ([], []) + results = self._filter_locations(results, genes, intervals, exclude_locations=exclude_locations, **kwargs) - if not (filter_qs or transcript_filters): - if any(val for key, val in (annotations or {}).items() if key != NEW_SV_FIELD) or any(val for val in (pathogenicity or {}).values()): - # Annotation filters restrict search to other dataset types - results = results.none() + if not (annotations_filter_q or transcript_filters): return results - filter_q = filter_qs[0] if filter_qs else None - for q in filter_qs[1:]: - filter_q |= q - if filter_q: - results = results.annotate(passes_annotation=filter_q) - filter_q = Q(passes_annotation=True) + if annotations_filter_q: + results = results.annotate(passes_annotation=annotations_filter_q) + annotations_filter_q = Q(passes_annotation=True) if transcript_filters: consequence_field = self.GENE_CONSEQUENCE_FIELD if genes else self.TRANSCRIPT_FIELD results = self._annotate_filtered_transcripts(results, consequence_field, transcript_filters, **kwargs) transcript_q = Q(filtered_transcript_consequences__not_empty=True) - if filter_q: - filter_q |= transcript_q + if annotations_filter_q: + annotations_filter_q |= transcript_q else: - filter_q = transcript_q + annotations_filter_q = transcript_q - return results.filter(filter_q) + return results.filter(annotations_filter_q) def _annotate_filtered_transcripts(self, results, consequence_field, transcript_filters, **kwargs): return results.annotate(**{ @@ -472,6 +469,16 @@ def _interval_query(self, chrom, start, end, offset=None, **kwargs): end = min(end + offset_pos, MAX_POS) return Q(xpos__range=(get_xpos(chrom, start), get_xpos(chrom, end))) + def get_parsed_annotations_filters(self, annotations=None, pathogenicity=None, **kwargs): + annotations = {k: v for k, v in (annotations or {}).items() if v} + pathogenicity = {k: v for k, v in (pathogenicity or {}).items() if v} + filter_qs, transcript_filters = self._parse_annotation_filters(annotations, pathogenicity) + filter_q = filter_qs[0] if filter_qs else None + for q in filter_qs[1:]: + filter_q |= q + + return {'annotations_filter_q': filter_q, 'transcript_filters': transcript_filters} + def _parse_annotation_filters(self, annotations, pathogenicity): raise NotImplementedError @@ -642,11 +649,7 @@ def _parse_annotation_filters(self, annotations, pathogenicity): filter_qs.append(Q(mitomapPathogenic=value)) elif field == SPLICE_AI_FIELD: if value and SPLICE_AI_FIELD in self.entry_model.PREDICTIONS: - pred_index = next ( - i for i, (pred_field, _) in enumerate(self.query.annotations['preds'].output_field.base_fields) - if pred_field == SPLICE_AI_FIELD - ) - filter_qs.append(Q(**{f'preds__{pred_index}__gte': float(value)})) + filter_qs.append(Q(preds__0__gte=float(value))) elif field not in SV_ANNOTATION_TYPES: allowed_consequences += value @@ -665,6 +668,9 @@ def _parse_annotation_filters(self, annotations, pathogenicity): if allowed_consequences: transcript_filters += self._allowed_consequences_filters(allowed_consequences) + if (annotations or pathogenicity) and not (filter_qs or transcript_filters): + raise InvalidDatasetTypeException + return filter_qs, transcript_filters def _allowed_consequences_filters(self, allowed_consequences): @@ -823,6 +829,10 @@ def _parse_annotation_filters(self, annotations, pathogenicity): if sv_type.startswith(self.model.SV_TYPE_FILTER_PREFIX) ])) + if not (filter_qs or transcript_filters) and (pathogenicity or any(key for key in annotations.keys() if key != NEW_SV_FIELD)): + # Annotation filters restrict search to other dataset types + raise InvalidDatasetTypeException + return filter_qs, transcript_filters def _filter_locations(self, results, genes, intervals, exclude_locations=False, padded_interval_end=None, **kwargs): @@ -1007,6 +1017,10 @@ def gt_stats_dict_rel(self): def genome_version(self): return self.annotations_model.ANNOTATION_CONSTANTS['genomeVersion'] + @property + def filtered_chrom(self): + return self.annotations_model.ANNOTATION_CONSTANTS.get('chrom') + @property def callset_filter_field(self): return 'callset' @@ -1016,7 +1030,7 @@ def _clinvar_path_q(cls, pathogenicity): return cls._clinvar_filter_q(pathogenicity, allowed_filters=CLINVAR_PATH_SIGNIFICANCES) def search(self, sample_data, exclude_keys=None, **kwargs): - entries = self.filter_locus(**kwargs) + entries = self if exclude_keys: entries = entries.exclude(key__in=exclude_keys) @@ -1072,9 +1086,6 @@ def _search_call_data(self, entries, sample_data, inheritance_mode=None, inherit elif exclude_projects: entries = entries.exclude(project_guid__in=exclude_projects) - if inheritance_mode in {X_LINKED_RECESSIVE, X_LINKED_RECESSIVE_MALE_AFFECTED}: - entries = entries.filter(self._interval_query('X', start=MIN_POS, end=MAX_POS)) - inheritance_q = None quality_q = None gt_filter = None @@ -1420,15 +1431,28 @@ def _family_passes_expression(pass_field): mapped_expression='x.1', output_field=models.ArrayField(models.StringField()), ) - def filter_locus(self, exclude_locations=False, require_gene_filter=False, intervals=None, genes=None, variant_ids=None, **kwargs): + def filter_locus(self, exclude_locations=False, require_gene_filter=False, intervals=None, genes=None, variant_ids=None, inheritance_mode=None, **kwargs): entries = self if variant_ids: entries = entries.filter_variant_ids(variant_ids) - if not (genes or intervals): - return entries + if inheritance_mode in {X_LINKED_RECESSIVE, X_LINKED_RECESSIVE_MALE_AFFECTED}: + if self.filtered_chrom and self.filtered_chrom != 'X': + raise InvalidDatasetTypeException + entries = entries.filter(self._interval_query('X', start=MIN_POS, end=MAX_POS)) + + if intervals and self.filtered_chrom: + intervals = [interval for interval in intervals if interval['chrom'] == self.filtered_chrom] + if not intervals and not exclude_locations: + raise InvalidDatasetTypeException + + if genes or intervals: + entries = self._filter_locations(entries, genes, intervals, exclude_locations=exclude_locations, require_gene_filter=require_gene_filter) + + return entries + def _filter_locations(self, entries, genes, intervals, exclude_locations=False, require_gene_filter=False): locus_q = None if genes: should_filter_interval = self._can_filter_gene_interval(genes) or exclude_locations @@ -1452,7 +1476,7 @@ def filter_locus(self, exclude_locations=False, require_gene_filter=False, inter return filter_func(locus_q) def _can_filter_gene_interval(self, genes): - return (not hasattr(self.model, 'geneId_ids')) or len(genes) < self.MAX_XPOS_FILTER_INTERVALS + return (not hasattr(self.model, 'geneId_ids')) or len(genes) < self.MAX_XPOS_FILTER_INTERVALS or self.filtered_chrom def search_padded_interval(self, chrom, pos, padding): interval_q = self._interval_query(chrom, start=max(pos - padding, MIN_POS), end=min(pos + padding, MAX_POS)) @@ -1618,6 +1642,8 @@ def filter_locus(self, *args, require_any_gene=False, parsed_variant_ids=None, i # improve performance by using the xpos projection intervals = [{'chrom': chrom, 'start': pos, 'end': pos} for chrom, pos, _, _ in parsed_variant_ids] variant_ids = ['-'.join([str(o) for o in variant_id]) for variant_id in parsed_variant_ids] + elif variant_ids: + raise InvalidDatasetTypeException entries = super().filter_locus(*args, intervals=intervals, genes=genes, variant_ids=variant_ids, **kwargs) @@ -1689,7 +1715,9 @@ def _annotate_calls(self, entries, *args, **kwargs): return entries - def filter_locus(self, *args, exclude_locations=False, intervals=None, genes=None, **kwargs): + def filter_locus(self, *args, exclude_locations=False, intervals=None, genes=None, exclude_svs=False, parsed_variant_ids=None, **kwargs): + if exclude_svs or parsed_variant_ids: + raise InvalidDatasetTypeException # SV interval filtering occurs after joining on annotations to correctly incorporate end position can_filter_gene_interval = self._can_filter_gene_interval(genes) if exclude_locations: diff --git a/clickhouse_search/models/search_models.py b/clickhouse_search/models/search_models.py index d3c609df9d..bcd47b2dcc 100644 --- a/clickhouse_search/models/search_models.py +++ b/clickhouse_search/models/search_models.py @@ -300,11 +300,11 @@ class BaseEntriesSnvIndel(BaseEntries): ('dp', models.UInt16Field(null=True, blank=True)), ] PREDICTIONS = { - 'dbnsfp': {}, - 'eigen': {}, 'splice_ai': { 'consequence_id': ('splice_ai_consequence', models.Enum8Field(choices=[(csq, csq_id) for csq_id, csq in BaseSpliceAi.CONSEQUENCE_CHOICES])), }, + 'dbnsfp': {}, + 'eigen': {}, } RANGE_PREDICTIONS = {} POPULATIONS = ['gnomad_exomes', 'gnomad_genomes', 'topmed'] diff --git a/clickhouse_search/search.py b/clickhouse_search/search.py index faa14f7a81..fd37ecb84b 100644 --- a/clickhouse_search/search.py +++ b/clickhouse_search/search.py @@ -10,6 +10,7 @@ from clickhouse_search.backend.fields import NamedTupleField from clickhouse_search.backend.functions import Array, ArrayFilter, ArrayIntersect, ArraySort, GroupArrayArray, If, Tuple, \ ArrayMap, Modulo +from clickhouse_search.managers import InvalidDatasetTypeException from clickhouse_search.models.gt_stats_models import PROJECT_GT_STATS_VIEW_CLASS_MAP from clickhouse_search.models.reference_data_models import BaseClinvar, BaseHgmd from clickhouse_search.models.search_models import BaseVariants, BaseVariantsSvGcnv, \ @@ -34,103 +35,189 @@ class InvalidSearchException(Exception): SELECTED_TRANSCRIPT_FIELD = 'selectedTranscript' -def get_clickhouse_variants(families, dataset_types, search, user, genome_version, sort=None, secondary_dataset_types=None, sample_data_by_dataset_type=None, encode_genotypes_json=False): +def get_clickhouse_variants(families, search, user, genome_version, sort=None, sample_data_by_dataset_type=None, encode_genotypes_json=False): inheritance_mode = search.get('inheritance_mode') has_comp_het = inheritance_mode in {RECESSIVE, COMPOUND_HET} has_x_chrom_comp_het = has_comp_het and _is_x_chrom_only(genome_version, **search) has_x_linked = inheritance_mode in {RECESSIVE, X_LINKED_RECESSIVE} and _has_x_chrom(genome_version, **search) - sample_data_by_dataset_type = sample_data_by_dataset_type or _get_sample_data( - families, - dataset_types, - secondary_dataset_types=secondary_dataset_types, - skip_multi_project_individual_guid=True, - annotate_affected_males=has_x_chrom_comp_het or has_x_linked, - has_comp_het=has_comp_het, - inheritance_mode=inheritance_mode, - inheritance_filter=search.get('inheritance_filter'), - allow_no_samples=bool(search.get('no_access_project_genome_version')), - ) + sample_data_by_dataset_type = sample_data_by_dataset_type or {} results = [] - family_guid = None exclude_keys = search.pop('exclude_keys', None) or {} exclude_key_pairs = search.pop('exclude_key_pairs', None) or {} - for dataset_type, sample_data in sample_data_by_dataset_type.items(): - logger.info(f'Loading {dataset_type} data for {sample_data["num_families"]} families', user) + searched_dataset_types = set() + sample_data_errors = set() + for dataset_type in ENTRY_CLASS_MAP[genome_version]: + try: + entry_qs, variants_qs, parsed_filters = _parse_dataset_type_query( + genome_version, dataset_type, families, sample_data_by_dataset_type, sample_data_errors, + annotate_affected_males=has_x_chrom_comp_het or has_x_linked, **search, + ) + except InvalidDatasetTypeException: + continue - family_guid = next(iter(next(iter(sample_data['sample_type_families'].values())))) + if dataset_type == Sample.DATASET_TYPE_VARIANT_CALLS and search.get('no_access_project_genome_version'): + logger.info('Looking up variants in projects with no user access', user) + results += _get_search_results( + entry_qs, variants_qs, **search, **parsed_filters, sample_data=None, skip_entry_fields=True, + exclude_projects=sample_data_by_dataset_type[dataset_type].get('project_guids'), + ) + searched_dataset_types.add(dataset_type) + + sample_data = sample_data_by_dataset_type[dataset_type] + if not sample_data: + continue + + logger.info(f'Loading {dataset_type} data for {sample_data["num_families"]} families', user) dataset_results = [] if inheritance_mode != COMPOUND_HET: - dataset_results += _get_search_results(genome_version, dataset_type, sample_data, exclude_keys=exclude_keys.get(dataset_type), **search) + dataset_results += _get_search_results( + entry_qs, variants_qs, sample_data, exclude_keys=exclude_keys.get(dataset_type), **search, **parsed_filters, + ) run_x_linked_male_search = has_x_linked and not (inheritance_mode == X_LINKED_RECESSIVE and sample_data.get('samples')) if run_x_linked_male_search: - affected_male_family_guids = { - s['family_guid'] for s in sample_data['samples'] if s['affected'] == AFFECTED and s['sex'] in MALE_SEXES - } if 'samples' in sample_data else sample_data['affected_male_family_guids'] - if affected_male_family_guids: - x_linked_sample_data = _affected_male_families(sample_data, affected_male_family_guids) - x_linked_search = {**search, 'inheritance_mode': X_LINKED_RECESSIVE_MALE_AFFECTED} - logger.info(f'Loading {dataset_type} X-linked male data for {x_linked_sample_data["num_families"]} families', user) - dataset_results += _get_search_results( - genome_version, dataset_type, x_linked_sample_data, exclude_keys=exclude_keys.get(dataset_type), - **x_linked_search, - ) + dataset_results += _get_x_linked_male_search_results( + entry_qs, variants_qs, dataset_type, user, sample_data, exclude_keys=exclude_keys.get(dataset_type), + **search, **parsed_filters, + ) if has_comp_het: comp_het_sample_data = sample_data if has_x_chrom_comp_het and 'affected_male_family_guids' in sample_data and dataset_type == Sample.DATASET_TYPE_VARIANT_CALLS: comp_het_sample_data = _no_affected_male_families(sample_data, user) - result_q = _get_data_type_comp_het_results_queryset(genome_version, dataset_type, comp_het_sample_data, exclude_key_pairs=exclude_key_pairs.get(dataset_type), **search) - dataset_results += _evaluate_results(result_q, is_comp_het=True) + result_q = _get_data_type_comp_het_results_queryset( + entry_qs, variants_qs, comp_het_sample_data, exclude_key_pairs=exclude_key_pairs.get(dataset_type), + **search, **parsed_filters, + ) + if result_q is not None: + dataset_results += _evaluate_results(result_q, is_comp_het=True) if 'samples' not in sample_data: _add_individual_guids(dataset_results, encode_genotypes_json=encode_genotypes_json) results += dataset_results + searched_dataset_types.add(dataset_type) - if has_comp_het and Sample.DATASET_TYPE_VARIANT_CALLS in sample_data_by_dataset_type and any( - dataset_type.startswith(Sample.DATASET_TYPE_SV_CALLS) for dataset_type in sample_data_by_dataset_type - ): - results += _get_multi_data_type_comp_het_results(genome_version, sample_data_by_dataset_type, user, exclude_key_pairs, **search) + if has_comp_het and any(dt.startswith(Sample.DATASET_TYPE_SV_CALLS) for dt in VARIANTS_CLASS_MAP[genome_version]): + results += _get_multi_data_type_comp_het_results(genome_version, families, sample_data_by_dataset_type, user, exclude_key_pairs, searched_dataset_types, **search) - if search.get('no_access_project_genome_version'): - logger.info('Looking up variants in projects with no user access', user) - results += _get_search_results( - genome_version, Sample.DATASET_TYPE_VARIANT_CALLS, **search, sample_data=None, skip_entry_fields=True, - exclude_projects=sample_data_by_dataset_type.get(Sample.DATASET_TYPE_VARIANT_CALLS, {}).get('project_guids'), - ) + if not searched_dataset_types: + _raise_dataset_type_errors(sample_data_errors, sample_data_by_dataset_type) logger.info(f'Total results: {len(results)}', user) - return get_sorted_search_results(results, sort, family_guid) + return get_sorted_search_results(results, sort, families) -def _get_search_results(genome_version, dataset_type, sample_data, skip_entry_fields=False, **search_kwargs): - entry_cls = ENTRY_CLASS_MAP[genome_version][dataset_type] - variants_cls = VARIANTS_CLASS_MAP[genome_version][dataset_type] - entries = entry_cls.objects.search(sample_data, **search_kwargs, skip_entry_fields=skip_entry_fields) - results = variants_cls.objects.subquery_join(entries).search(**search_kwargs, skip_entry_fields=skip_entry_fields) +def _parse_dataset_type_query(genome_version, dataset_type, families, sample_data_by_dataset_type, sample_data_errors, annotate_affected_males, inheritance_mode=None, inheritance_filter=None, no_access_project_genome_version=None, **search_kwargs): + entry_qs = ENTRY_CLASS_MAP[genome_version][dataset_type].objects.filter_locus(inheritance_mode=inheritance_mode, **search_kwargs) + variants_qs = VARIANTS_CLASS_MAP[genome_version][dataset_type].objects + parsed_filters = variants_qs.get_parsed_annotations_filters(**search_kwargs) + + if dataset_type in sample_data_by_dataset_type: + sample_data = sample_data_by_dataset_type[dataset_type] + else: + sample_data = _get_sample_data( + families, + dataset_type, + inheritance_mode=inheritance_mode, + inheritance_filter=inheritance_filter, + annotate_affected_males=annotate_affected_males, + allow_no_samples=bool(no_access_project_genome_version), + ) + sample_data_by_dataset_type[dataset_type] = sample_data + + individual_affected_status = (inheritance_filter or {}).get('affected') + genotype_filter = (inheritance_filter or {}).get('genotype') + if sample_data and inheritance_mode and (not sample_data['num_families'] or (individual_affected_status and all( + (individual_affected_status.get(s['individual_guid']) or s['affected']) != AFFECTED for s in sample_data.get('samples', []) + ))): + sample_data_errors.add( + 'Inheritance based search is disabled in families with no data loaded for affected individuals') + sample_data_by_dataset_type[dataset_type] = None + if sample_data and genotype_filter and all(s['individual_guid'] not in genotype_filter for s in sample_data.get('samples', [])): + sample_data_errors.add('Invalid custom inheritance') + sample_data_by_dataset_type[dataset_type] = None + + return entry_qs, variants_qs, parsed_filters + + +def _raise_dataset_type_errors(sample_data_errors, sample_data_by_dataset_type): + if sample_data_errors: + raise InvalidSearchException(next(iter(sample_data_errors))) + no_data_type = next(data_type for data_type, data in sample_data_by_dataset_type.items() if not data) + if no_data_type.startswith(Sample.DATASET_TYPE_SV_CALLS): + no_data_type = Sample.DATASET_TYPE_SV_CALLS + raise InvalidSearchException(f'Unable to search against dataset type "{no_data_type}"') + + +def _get_search_results(entry_qs, variants_qs, sample_data, skip_entry_fields=False, **search_kwargs): + entries = entry_qs.search(sample_data, skip_entry_fields=skip_entry_fields, **search_kwargs) + results = variants_qs.subquery_join(entries).search(skip_entry_fields=skip_entry_fields, **search_kwargs) return _evaluate_results(results.result_values(skip_entry_fields=skip_entry_fields)) +def _get_x_linked_male_search_results(entry_qs, variants_qs, dataset_type, user, sample_data, inheritance_mode=None, **search_kwargs): + affected_male_family_guids = { + s['family_guid'] for s in sample_data['samples'] if s['affected'] == AFFECTED and s['sex'] in MALE_SEXES + } if 'samples' in sample_data else sample_data['affected_male_family_guids'] + if not affected_male_family_guids: + return [] + try: + x_linked_entry_qs = entry_qs.filter_locus(inheritance_mode=X_LINKED_RECESSIVE_MALE_AFFECTED) + except InvalidDatasetTypeException: + return [] + + x_linked_sample_data = _affected_male_families(sample_data, affected_male_family_guids) + logger.info(f'Loading {dataset_type} X-linked male data for {x_linked_sample_data["num_families"]} families', user) + return _get_search_results( + x_linked_entry_qs, variants_qs, x_linked_sample_data, inheritance_mode=X_LINKED_RECESSIVE_MALE_AFFECTED, **search_kwargs, + ) + + def _evaluate_results(result_q, is_comp_het=False): results = [list(result[1:]) if is_comp_het else result for result in result_q[:MAX_VARIANTS + 1]] if len(results) > MAX_VARIANTS: raise InvalidSearchException('This search returned too many results') return results -def _get_multi_data_type_comp_het_results(genome_version, sample_data_by_dataset_type, user, exclude_key_pairs, annotations=None, annotations_secondary=None, inheritance_mode=None, **search_kwargs): +def _get_multi_data_type_comp_het_results(genome_version, all_families, sample_data_by_dataset_type, user, exclude_key_pairs, searched_dataset_types, annotations=None, annotations_secondary=None, inheritance_mode=None, **search_kwargs): if annotations_secondary: annotations = { **annotations, **{k: v + annotations[k] if k in annotations else v for k, v in annotations_secondary.items()}, } + try: + snv_indel_entry_qs = ENTRY_CLASS_MAP[genome_version][Sample.DATASET_TYPE_VARIANT_CALLS].objects.filter_locus(**search_kwargs) + snv_indel_variants_qs = VARIANTS_CLASS_MAP[genome_version][Sample.DATASET_TYPE_VARIANT_CALLS].objects + snv_indel_parsed_filters = snv_indel_variants_qs.get_parsed_annotations_filters(annotations=annotations, **search_kwargs) + sv_variants_cls = VARIANTS_CLASS_MAP[genome_version][f'{Sample.DATASET_TYPE_SV_CALLS}_{Sample.SAMPLE_TYPE_WES}'] + sv_parsed_filters = sv_variants_cls.objects.get_parsed_annotations_filters(annotations=annotations, **search_kwargs) + except InvalidDatasetTypeException: + return [] + + if Sample.DATASET_TYPE_VARIANT_CALLS not in sample_data_by_dataset_type: + sample_data_by_dataset_type[Sample.DATASET_TYPE_VARIANT_CALLS] = _get_sample_data( + all_families, + Sample.DATASET_TYPE_VARIANT_CALLS, + inheritance_mode=COMPOUND_HET, + inheritance_filter=search_kwargs.get('inheritance_filter'), + ) snv_indel_sample_data = sample_data_by_dataset_type[Sample.DATASET_TYPE_VARIANT_CALLS] + if not (snv_indel_sample_data or {}).get('num_families'): + return [] snv_indel_families = set().union(*snv_indel_sample_data['sample_type_families'].values()) results = [] for sample_type in [Sample.SAMPLE_TYPE_WES, Sample.SAMPLE_TYPE_WGS]: sv_dataset_type = f'{Sample.DATASET_TYPE_SV_CALLS}_{sample_type}' - sv_sample_data = sample_data_by_dataset_type.get(sv_dataset_type, {}) + if sv_dataset_type not in sample_data_by_dataset_type: + sample_data_by_dataset_type[sv_dataset_type] = _get_sample_data( + all_families, + sv_dataset_type, + inheritance_mode=COMPOUND_HET, + inheritance_filter=search_kwargs.get('inheritance_filter'), + ) + sv_sample_data = sample_data_by_dataset_type[sv_dataset_type] or {} sv_families = set(sv_sample_data.get('sample_type_families', {}).get(sample_type, [])) families = snv_indel_families.intersection(sv_families) if not families: @@ -156,55 +243,57 @@ def _get_multi_data_type_comp_het_results(genome_version, sample_data_by_dataset 'samples': [s for s in sv_sample_data['samples'] if s['family_guid'] in families] if 'samples' in sv_sample_data else None, } - entries = ENTRY_CLASS_MAP[genome_version][Sample.DATASET_TYPE_VARIANT_CALLS].objects.search( + entries = snv_indel_entry_qs.search( type_snv_indel_sample_data, inheritance_mode=COMPOUND_HET_ALLOW_HOM_ALTS, annotate_carriers=True, - annotate_hom_alts=True, annotations=annotations, **search_kwargs, + annotate_hom_alts=True, **search_kwargs, ) - variants_cls = VARIANTS_CLASS_MAP[genome_version][Sample.DATASET_TYPE_VARIANT_CALLS] - snv_indel_q = variants_cls.objects.subquery_join(entries).search(annotations=annotations, **search_kwargs) + snv_indel_q = snv_indel_variants_qs.subquery_join(entries).search(**search_kwargs, **snv_indel_parsed_filters) - sv_entries = ENTRY_CLASS_MAP[genome_version][sv_dataset_type].objects.search( - sv_sample_data, annotations=annotations, **search_kwargs, inheritance_mode=COMPOUND_HET, annotate_carriers=True, + sv_entries = ENTRY_CLASS_MAP[genome_version][sv_dataset_type].objects.filter_locus(**search_kwargs).search( + sv_sample_data, **search_kwargs, inheritance_mode=COMPOUND_HET, annotate_carriers=True, ) sv_variants_cls = VARIANTS_CLASS_MAP[genome_version][sv_dataset_type] - sv_q = sv_variants_cls.objects.subquery_join(sv_entries).search(annotations=annotations, **search_kwargs) + sv_q = sv_variants_cls.objects.subquery_join(sv_entries).search(**search_kwargs, **sv_parsed_filters) result_q = _get_comp_het_results_queryset( - variants_cls, snv_indel_q, sv_q, len(families), exclude_key_pairs=exclude_key_pairs.get(f'{Sample.DATASET_TYPE_VARIANT_CALLS},{sv_dataset_type}'), + snv_indel_variants_qs, snv_indel_q, sv_q, len(families), + exclude_key_pairs=exclude_key_pairs.get(f'{Sample.DATASET_TYPE_VARIANT_CALLS},{sv_dataset_type}'), ) dataset_results = _evaluate_results(result_q, is_comp_het=True) if not sv_sample_data['samples']: _add_individual_guids(dataset_results) results += dataset_results + searched_dataset_types.add(sv_dataset_type) + + if not searched_dataset_types: + raise InvalidSearchException( + 'Unable to search for comp-het pairs with dataset type "SV". This may be because inheritance based search is disabled in families with no loaded affected individuals' + ) return results -def _get_data_type_comp_het_results_queryset(genome_version, dataset_type, sample_data, annotations=None, annotations_secondary=None, pathogenicity=None, inheritance_mode=None, exclude_key_pairs=None, split_pathogenicity_annotations=False, **search_kwargs): - entry_cls = ENTRY_CLASS_MAP[genome_version][dataset_type] - variants_cls = VARIANTS_CLASS_MAP[genome_version][dataset_type] - entries = entry_cls.objects.search( - sample_data, **search_kwargs, inheritance_mode=COMPOUND_HET, pathogenicity=pathogenicity, annotations=annotations, annotate_carriers=True, - ) +def _get_data_type_comp_het_results_queryset(entry_qs, variants_qs, sample_data, annotations_secondary=None, inheritance_mode=None, exclude_key_pairs=None, **search_kwargs): + parsed_secondary_filters = {} + if annotations_secondary: + try: + parsed_secondary_filters = variants_qs.get_parsed_annotations_filters(annotations=annotations_secondary) + except InvalidDatasetTypeException: + return None - if split_pathogenicity_annotations: - pathogenicity_secondary = pathogenicity - pathogenicity = None - else: - annotations_secondary = annotations_secondary or annotations - pathogenicity_secondary = pathogenicity - primary_q = variants_cls.objects.subquery_join(entries).search( - annotations=annotations, pathogenicity=pathogenicity, **search_kwargs, - ) - secondary_q = variants_cls.objects.subquery_join(entries).search( - annotations=annotations_secondary, pathogenicity=pathogenicity_secondary, **search_kwargs, + entries = entry_qs.search( + sample_data, **search_kwargs, inheritance_mode=COMPOUND_HET, annotate_carriers=True, ) - return _get_comp_het_results_queryset(variants_cls, primary_q, secondary_q, sample_data['num_families'], exclude_key_pairs) + primary_q = variants_qs.subquery_join(entries).search(**search_kwargs) + secondary_kwargs = {**search_kwargs, **parsed_secondary_filters} + secondary_q = variants_qs.subquery_join(entries).search(**secondary_kwargs) + + return _get_comp_het_results_queryset(variants_qs, primary_q, secondary_q, sample_data['num_families'], exclude_key_pairs) -def _get_comp_het_results_queryset(variants_cls, primary_q, secondary_q, num_families, exclude_key_pairs): - results = variants_cls.objects.search_compound_hets(primary_q, secondary_q) +def _get_comp_het_results_queryset(variants_qs, primary_q, secondary_q, num_families, exclude_key_pairs): + results = variants_qs.search_compound_hets(primary_q, secondary_q) if results.has_annotation('primary_carriers') and results.has_annotation('secondary_carriers'): results = results.annotate( @@ -316,8 +405,8 @@ def _set_individual_guids(result, sample_map, encode_genotypes_json): result['genotypes'] = genotypes -def get_sorted_search_results(results, sort, family_guid): - sort_metadata = _get_sort_gene_metadata(sort, results, family_guid) +def get_sorted_search_results(results, sort, families): + sort_metadata = _get_sort_gene_metadata(sort, results, families) sort_key = _get_sort_key(sort, sort_metadata) sorted_results = sorted([ sorted(result, key=sort_key) if isinstance(result, list) else result for result in results @@ -412,21 +501,18 @@ def _is_matched_minimal_transcript(transcript, minimal_transcript): and transcript.get('spliceregion', {}).get('extended_intronic_splice_region_variant') == minimal_transcript.get('extendedIntronicSpliceRegionVariant')) -def _get_valid_samples(families, dataset_types, secondary_dataset_types, allow_no_samples): +def _get_valid_samples(families, dataset_type, sample_type, allow_no_samples): samples = Sample.objects.filter(individual__family__in=families, is_active=True) if not samples.exists(): if allow_no_samples: return None raise InvalidSearchException(f'No search data found for families {", ".join([f.family_id for f in families])}') - if dataset_types: - samples = samples.filter(dataset_type__in={*dataset_types, *(secondary_dataset_types or [])}) - if not samples.exists(): - if allow_no_samples: - return None - raise InvalidSearchException(f'Unable to search against dataset type "{dataset_types[0]}"') + samples = samples.filter(dataset_type=dataset_type) + if sample_type: + samples = samples.filter(sample_type=sample_type) - mismatch_affected_samples = samples.values('sample_id', 'dataset_type').annotate( + mismatch_affected_samples = samples.values('sample_id').annotate( projects=ArrayAgg('individual__family__project__name', distinct=True), affected=ArrayAgg('individual__affected', distinct=True), ).filter(affected__len__gt=1) @@ -439,10 +525,8 @@ def _get_valid_samples(families, dataset_types, secondary_dataset_types, allow_n return samples -def _get_grouped_samples(samples, skip_multi_project_individual_guid, affected_family_only, annotate_affected_males): - skip_individual_guid = ( - skip_multi_project_individual_guid and samples.values('individual__family__project_id').distinct().count() > 1 - ) +def _get_sample_metadata(samples, affected_family_only, annotate_affected_males): + skip_individual_guid = samples.values('individual__family__project_id').distinct().count() > 1 family_array_kwargs = {'distinct': True} if affected_family_only: family_array_kwargs['filter'] = Q(individual__affected=Individual.AFFECTED_STATUS_AFFECTED) @@ -464,108 +548,70 @@ def _get_grouped_samples(samples, skip_multi_project_individual_guid, affected_f family_guid=F('individual__family__guid'), individual_guid=F('individual__guid'), )) - sample_data = samples.values('dataset_type', 'sample_type').annotate(**annotations) - if affected_family_only: - sample_data = sample_data.filter(family_guids__len__gt=0) - return sample_data + return samples.aggregate(**annotations) -def _get_sample_data(families, dataset_types, skip_multi_project_individual_guid=False, annotate_affected_males=False, allow_no_samples=False, has_comp_het=False, secondary_dataset_types=None, inheritance_mode=None, inheritance_filter=None): - samples = _get_valid_samples(families, dataset_types, secondary_dataset_types, allow_no_samples) +def _get_sample_data(families, dataset_type, annotate_affected_males=False, allow_no_samples=False, inheritance_mode=None, inheritance_filter=None): + sample_type = None + if dataset_type.startswith(Sample.DATASET_TYPE_SV_CALLS): + dataset_type, sample_type = dataset_type.split('_') + samples = _get_valid_samples(families, dataset_type, sample_type, allow_no_samples) if not samples: return {} individual_affected_status = (inheritance_filter or {}).get('affected') - genotype_filter = (inheritance_filter or {}).get('genotype') affected_family_only = inheritance_mode and not individual_affected_status - sample_data = _get_grouped_samples(samples, skip_multi_project_individual_guid, affected_family_only, annotate_affected_males) - - samples_by_dataset_type = {} - for data in sample_data: - dataset_type = data.pop('dataset_type') - sample_type = data.pop('sample_type') - if dataset_type == Sample.DATASET_TYPE_SV_CALLS: - dataset_type = f'{Sample.DATASET_TYPE_SV_CALLS}_{sample_type}' + sample_data = _get_sample_metadata(samples, affected_family_only, annotate_affected_males) - if genotype_filter and all(s['individual_guid'] not in genotype_filter for s in data.get('samples', [])): - continue - if inheritance_mode and individual_affected_status and all( - (individual_affected_status.get(s['individual_guid']) or s['affected']) != AFFECTED for s in data.get('samples', []) - ): - continue - - family_guids = set(data.pop('family_guids')) - if dataset_type in samples_by_dataset_type: - other_type_data = samples_by_dataset_type[dataset_type] - other_sample_type, other_type_family_guids = next(iter(other_type_data['sample_type_families'].items())) - sample_type_families = { - other_sample_type: other_type_family_guids - family_guids, - sample_type: family_guids - other_type_family_guids, - 'multi': family_guids.intersection(other_type_family_guids), - } - if sample_type_families['multi']: - data['family_missing_type_samples'] = defaultdict(lambda: defaultdict(list)) - individual_samples = {s['individual_guid']: s for s in data.get('samples', [])} - other_type_individual_samples = {s['individual_guid']: s for s in other_type_data.get('samples', [])} - for individual_guid in set(individual_samples.keys()) - set(other_type_individual_samples.keys()): - sample = individual_samples[individual_guid] - data['family_missing_type_samples'][sample['family_guid']][other_sample_type].append(sample['sample_id']) - for individual_guid in set(other_type_individual_samples.keys()) - set(individual_samples.keys()): - sample = other_type_individual_samples[individual_guid] - data['family_missing_type_samples'][sample['family_guid']][sample_type].append(sample['sample_id']) - data['sample_type_families'] = {k: v for k, v in sample_type_families.items() if v} - for key in ['project_guids', 'samples']: - if key in data: - data[key] += other_type_data[key] + family_guids = set(sample_data.pop('family_guids')) + sample_data['num_families'] = len(family_guids) + if sample_type: + sample_data['sample_type_families'] = {sample_type: family_guids} + else: + if sample_data.get('samples'): + sample_data['sample_type_families'] = defaultdict(set) + for sample in sample_data['samples']: + sample_data['sample_type_families'][sample['sample_type']].add(sample['family_guid']) else: - data['sample_type_families'] = {sample_type: family_guids} - - samples_by_dataset_type[dataset_type] = { - **data, - 'num_families': len(set().union(*data['sample_type_families'].values())), - } - - _validate_dataset_types( - samples_by_dataset_type, dataset_types, secondary_dataset_types, has_comp_het, genotype_filter, - ) - _add_missing_multi_type_samples(samples, samples_by_dataset_type) - - return samples_by_dataset_type - - -def _validate_dataset_types(samples_by_dataset_type, dataset_types, secondary_dataset_types, has_comp_het, genotype_filter): - if not samples_by_dataset_type: - raise InvalidSearchException( - 'Invalid custom inheritance' if genotype_filter else - 'Inheritance based search is disabled in families with no data loaded for affected individuals' - ) - - if has_comp_het and secondary_dataset_types: - sample_dataset_types = { - Sample.DATASET_TYPE_SV_CALLS if dt.startswith(Sample.DATASET_TYPE_SV_CALLS) else dt - for dt in samples_by_dataset_type - } - invalid_type = next(( - dts[0] for dts in [dataset_types, secondary_dataset_types] if dts and not sample_dataset_types.intersection(dts) - ), None) - if invalid_type: - raise InvalidSearchException( - f'Unable to search for comp-het pairs with dataset type "{invalid_type}". This may be because inheritance based search is disabled in families with no loaded affected individuals' + sample_data['sample_type_families'] = dict( + samples.filter(individual__family__guid__in=family_guids).values('sample_type').values_list( + 'sample_type', ArrayAgg('individual__family__guid', distinct=True), + ) ) + if len(sample_data['sample_type_families']) == 2: + st_families_1, st_families_2 = sample_data['sample_type_families'].values() + multi_families = set(st_families_1).intersection(st_families_2) + if multi_families: + sample_data['sample_type_families'] = { + sample_type: set(families) - multi_families + for sample_type, families in sample_data['sample_type_families'].items() + if set(families) - multi_families + } + sample_data['sample_type_families']['multi'] = multi_families + _add_missing_multi_type_samples(samples, sample_data) + + return sample_data -def _add_missing_multi_type_samples(samples, samples_by_dataset_type): - for dataset_type, data in samples_by_dataset_type.items(): - if data['sample_type_families'].get('multi') and 'samples' not in data: - individual_samples = samples.filter( - dataset_type=dataset_type, individual__family__guid__in=data['sample_type_families']['multi'], - ).values('individual_id', family_guid=F('individual__family__guid')).annotate( - samples=ArrayAgg(JSONObject(sample_id='sample_id', sample_type='sample_type')) - ).filter(samples__len=1) - for agg in individual_samples: - sample = agg['samples'][0] - missing_type = Sample.SAMPLE_TYPE_WES if sample['sample_type'] == Sample.SAMPLE_TYPE_WGS else Sample.SAMPLE_TYPE_WGS - data['family_missing_type_samples'][agg['family_guid']][missing_type].append(sample['sample_id']) +def _add_missing_multi_type_samples(samples, data): + data['family_missing_type_samples'] = defaultdict(lambda: defaultdict(list)) + if 'samples' in data: + individual_sample_types = defaultdict(list) + for s in data['samples']: + individual_sample_types[s['individual_guid']].append(s) + individual_samples = [ + {'samples': samples, 'family_guid': samples[0]['family_guid']} + for samples in individual_sample_types.values() if len(samples) == 1 + ] + else: + individual_samples = samples.filter(individual__family__guid__in=data['sample_type_families']['multi'], + ).values('individual_id', family_guid=F('individual__family__guid')).annotate( + samples=ArrayAgg(JSONObject(sample_id='sample_id', sample_type='sample_type')) + ).filter(samples__len=1) + for agg in individual_samples: + sample = agg['samples'][0] + missing_type = Sample.SAMPLE_TYPE_WES if sample['sample_type'] == Sample.SAMPLE_TYPE_WGS else Sample.SAMPLE_TYPE_WGS + data['family_missing_type_samples'][agg['family_guid']][missing_type].append(sample['sample_id']) def _no_affected_male_families(sample_data, user): @@ -617,14 +663,14 @@ def _has_x_chrom(genome_version, genes=None, intervals=None, **kwargs): OMIM_SORT: lambda gene_ids, _: set(Omim.objects.filter( gene__gene_id__in=gene_ids, phenotype_mim_number__isnull=False, ).values_list('gene__gene_id', flat=True)), - PRIORITIZED_GENE_SORT: lambda gene_ids, family_guid: { + PRIORITIZED_GENE_SORT: lambda gene_ids, families: { agg['gene_id']: agg['min_rank'] for agg in PhenotypePrioritization.objects.filter( - gene_id__in=gene_ids, individual__family__guid=family_guid, rank__lte=100, + gene_id__in=gene_ids, individual__family_id=families[0].id, rank__lte=100, ).values('gene_id').annotate(min_rank=Min('rank')) }, } -def _get_sort_gene_metadata(sort, results, family_guid): +def _get_sort_gene_metadata(sort, results, families): get_metadata = GENE_SORTS.get(sort) if not get_metadata: return None @@ -638,7 +684,7 @@ def _get_sort_gene_metadata(sort, results, family_guid): gene_ids.update([t['geneId'] for t in variant[TRANSCRIPT_CONSEQUENCES_FIELD]]) else: gene_ids.update(variant.get('transcripts', {}).keys()) - return get_metadata(gene_ids, family_guid) + return get_metadata(gene_ids, families) MAX_SORT_RANK = 1e10 @@ -742,13 +788,9 @@ def _get_sort_key(sort, gene_metadata): return lambda x: tuple(expr(x[0] if isinstance(x, list) else x) for expr in [*sort_expressions, lambda x: x[XPOS_SORT_KEY]]) -def _clickhouse_variant_lookup(variant_id, parsed_variant_id, genome_version, data_type, affected_only=False, hom_only=False): - entry_cls = ENTRY_CLASS_MAP[genome_version][data_type] +def _clickhouse_variant_lookup(entries, genome_version, data_type, affected_only=False, hom_only=False): variants_cls = VARIANTS_CLASS_MAP[genome_version][data_type] - entries = entry_cls.objects.filter_locus( - variant_ids=[variant_id], parsed_variant_ids=[parsed_variant_id] if parsed_variant_id else [], - ) entries = _filter_lookup_entries(entries, affected_only, hom_only) entries = entries.result_values() results = variants_cls.objects.subquery_join(entries) @@ -767,25 +809,43 @@ def _filter_lookup_entries(entries, affected_only, hom_only): entries = entries.filter(calls__array_exists={'gt': (2,)}) return entries -def clickhouse_variant_lookup(user, variant_id, parsed_variant_id, dataset_type, sample_type, genome_version, affected_only, hom_only): - is_sv = dataset_type == Sample.DATASET_TYPE_SV_CALLS - data_type = f'{dataset_type}_{sample_type}' if is_sv else dataset_type +def clickhouse_variant_lookup(user, variant_id, parsed_variant_id, sample_type, genome_version, affected_only, hom_only): + data_type = entry_qs = None + for dataset_type, entry_cls in sorted(ENTRY_CLASS_MAP[genome_version].items()): + try: + entry_qs = entry_cls.objects.filter_locus( + variant_ids=[variant_id], parsed_variant_ids=[parsed_variant_id] if parsed_variant_id else [], + ) + except InvalidDatasetTypeException: + continue + if dataset_type.startswith(Sample.DATASET_TYPE_SV_CALLS): + if not sample_type: + raise InvalidSearchException('Sample type must be specified to look up a structural variant') + elif not dataset_type.endswith(sample_type): + continue + data_type = dataset_type + break + if data_type is None: + raise InvalidSearchException('Invalid genome build for dataset type') + logger.info(f'Looking up variant {variant_id} with data type {data_type}', user) variant = _clickhouse_variant_lookup( - variant_id, parsed_variant_id, genome_version, data_type, affected_only=affected_only, hom_only=hom_only, + entry_qs, genome_version, data_type, affected_only=affected_only, hom_only=hom_only, ) if variant: _add_liftover_genotypes(variant, data_type, parsed_variant_id, affected_only, hom_only) else: lifted_genome_version = next(gv for gv in ENTRY_CLASS_MAP.keys() if gv != genome_version) - if ENTRY_CLASS_MAP[lifted_genome_version].get(data_type): + lifted_entry_cls = ENTRY_CLASS_MAP[lifted_genome_version].get(data_type) + if lifted_entry_cls: from seqr.utils.search.utils import run_liftover liftover_results = run_liftover(lifted_genome_version, parsed_variant_id[0], parsed_variant_id[1]) if liftover_results: lifted_id = (liftover_results[0], liftover_results[1], *parsed_variant_id[2:]) + entry_qs = lifted_entry_cls.objects.filter_locus(parsed_variant_ids=[lifted_id]) variant = _clickhouse_variant_lookup( - None, lifted_id, lifted_genome_version, data_type, affected_only=affected_only, hom_only=hom_only, + entry_qs, lifted_genome_version, data_type, affected_only=affected_only, hom_only=hom_only, ) if not variant: @@ -793,7 +853,7 @@ def clickhouse_variant_lookup(user, variant_id, parsed_variant_id, dataset_type, variants = [variant] - if is_sv and variant['svType'] in {'DEL', 'DUP'}: + if variant.get('svType') in {'DEL', 'DUP'}: other_sample_type, other_entry_class = next( (dt, cls) for dt, cls in ENTRY_CLASS_MAP[genome_version].items() if dt != data_type and dt.startswith(Sample.DATASET_TYPE_SV_CALLS) @@ -804,8 +864,9 @@ def clickhouse_variant_lookup(user, variant_id, parsed_variant_id, dataset_type, entries = other_entry_class.objects.search_padded_interval(variant['chrom'], variant['pos'], padding) entries = _filter_lookup_entries(entries, affected_only, hom_only) results = other_variants_cls.objects.subquery_join(entries).search( - padded_interval_end=(variant['end'], padding), - annotations={'structural': [variant['svType'], f"gCNV_{variant['svType']}"]}, + padded_interval_end=(variant['end'], padding), **other_variants_cls.objects.get_parsed_annotations_filters( + annotations={'structural': [variant['svType'], f"gCNV_{variant['svType']}"]}, + ), ) variants += list(results.result_values()) @@ -827,7 +888,7 @@ def _add_liftover_genotypes(variant, data_type, variant_id, affected_only, hom_o def get_clickhouse_genotypes(project_guid, family_guids, genome_version, dataset_type, keys, additional_fields=None): - sample_data = _get_sample_data(Family.objects.filter(guid__in=family_guids), None)[dataset_type] + sample_data = _get_sample_data(Family.objects.filter(guid__in=family_guids), dataset_type) entries = ENTRY_CLASS_MAP[genome_version][dataset_type].objects.filter( project_guid=project_guid, family_guid__in=family_guids, key__in=keys, ) diff --git a/clickhouse_search/search_tests.py b/clickhouse_search/search_tests.py index bd6ad73664..a5078a8996 100644 --- a/clickhouse_search/search_tests.py +++ b/clickhouse_search/search_tests.py @@ -1254,15 +1254,10 @@ def test_variant_lookup(self, mock_liftover): individual_guids=['I000004_hg00731', 'I000005_hg00732', 'I000006_hg00733'], ) - url = url.replace('38', '37') + url = url.replace('M-4429-G-A', 'phase2_DEL_chr14_4640').replace('38', '37') response = self.client.get(url) self.assertEqual(response.status_code, 400) - self.assertEqual(response.json(), {'error': 'MITO variants are not available for GRCh37'}) - - url = url.replace('M-4429-G-A', 'phase2_DEL_chr14_4640') - response = self.client.get(url) - self.assertEqual(response.status_code, 400) - self.assertEqual(response.json(), {'error': 'SV variants are not available for GRCh37'}) + self.assertEqual(response.json(), {'error': 'Invalid genome build for dataset type'}) url = url.replace('37', '38') response = self.client.get(url) diff --git a/seqr/management/commands/tag_seqr_prioritized_variants.py b/seqr/management/commands/tag_seqr_prioritized_variants.py index 8e4e63fc88..71b4a1cd9b 100644 --- a/seqr/management/commands/tag_seqr_prioritized_variants.py +++ b/seqr/management/commands/tag_seqr_prioritized_variants.py @@ -5,7 +5,7 @@ from django.db.models import Q from django.db.models.functions import JSONObject -from clickhouse_search.search import get_clickhouse_variants +from clickhouse_search.search import get_clickhouse_variants, ENTRY_CLASS_MAP from panelapp.models import PaLocusListGene from reference_data.models import GENOME_VERSION_GRCh38 from seqr.models import Project, Family, Individual, Sample, LocusList @@ -191,8 +191,8 @@ }, 'Clinvar Pathogenic - Compound Heterozygous': { 'inheritance_mode': COMPOUND_HET, - 'split_pathogenicity_annotations': True, - 'annotations': HIGH_MODERATE_ANNOTATIONS, + 'annotations': {}, + 'annotations_secondary': HIGH_MODERATE_ANNOTATIONS, **CLINVAR_RECESSIVE_SEARCH, }, 'Clinvar Both Pathogenic - Compound Heterozygous': { @@ -236,8 +236,8 @@ }, 'Compound Heterozygous - Clinvar Pathogenic/ High Splice AI': { 'inheritance_mode': COMPOUND_HET, - 'split_pathogenicity_annotations': True, - 'annotations': { + 'annotations': {}, + 'annotations_secondary': { 'splice_ai': 0.5, }, **CLINVAR_RECESSIVE_SEARCH, @@ -628,8 +628,10 @@ def _run_multi_data_type_comp_het_search(cls, updates, search_counts, samples_by @staticmethod def _execute_search(sample_data_by_dataset_type, search_name, family_variant_data, family_guid_map, **kwargs): results = get_clickhouse_variants( - families=None, dataset_types=None, search=kwargs, user=None, genome_version=GENOME_VERSION_GRCh38, - sample_data_by_dataset_type=sample_data_by_dataset_type, encode_genotypes_json=True + families=None, search=kwargs, user=None, genome_version=GENOME_VERSION_GRCh38, + encode_genotypes_json=True, sample_data_by_dataset_type={ + **{dt: None for dt in ENTRY_CLASS_MAP[GENOME_VERSION_GRCh38]}, **sample_data_by_dataset_type, + }, ) for result in results: if isinstance(result, list): diff --git a/seqr/management/tests/tag_seqr_prioritized_variants_tests.py b/seqr/management/tests/tag_seqr_prioritized_variants_tests.py index 6519a3bd25..f03bd95a59 100644 --- a/seqr/management/tests/tag_seqr_prioritized_variants_tests.py +++ b/seqr/management/tests/tag_seqr_prioritized_variants_tests.py @@ -151,7 +151,7 @@ def _assert_expected_logs(self, num_new=0, num_unchanged=0, creation_stats=None) 'MITO', 1, MITO_MATCHES, **creation_stats.get('MITO', {}), ) + self._dataset_type_logs( 'multi data type', 1, MULTI_TYPE_MATCHES, **creation_stats.get('MULTI', {}), - search_dataset_types=['SNV_INDEL', 'SV_WES', 'SNV_INDEL/SV_WES'], + search_dataset_types=['SNV_INDEL', 'SNV_INDEL/SV_WES'], ) + [ (f'Tagged {num_new} new and 0 previously tagged variants in 1 families, found {num_unchanged} unchanged tags:', None) ] + [(f' {criteria}: {count} variants', None) for criteria, (count, _) in SNV_INDEL_MATCHES.items()] + [ diff --git a/seqr/migrations/0086_alter_rnasample_tissue_type.py b/seqr/migrations/0086_alter_rnasample_tissue_type.py new file mode 100644 index 0000000000..09da7b827f --- /dev/null +++ b/seqr/migrations/0086_alter_rnasample_tissue_type.py @@ -0,0 +1,18 @@ +# Generated by Django 4.2.24 on 2026-03-30 19:17 + +from django.db import migrations, models + + +class Migration(migrations.Migration): + + dependencies = [ + ('seqr', '0085_savedvariant_main_transcript_savedvariant_sv_type'), + ] + + operations = [ + migrations.AlterField( + model_name='rnasample', + name='tissue_type', + field=models.CharField(choices=[('WB', 'whole_blood'), ('F', 'fibroblasts'), ('M', 'muscle'), ('L', 'lymphocytes'), ('A', 'airway_cultured_epithelium'), ('B', 'brain'), ('LV', 'liver'), ('S', 'skin')], max_length=2), + ), + ] diff --git a/seqr/models.py b/seqr/models.py index 36363d728c..e7ab32f379 100644 --- a/seqr/models.py +++ b/seqr/models.py @@ -767,6 +767,8 @@ class RnaSample(ModelWithGUID): ('L', 'lymphocytes'), ('A', 'airway_cultured_epithelium'), ('B', 'brain'), + ('LV', 'liver'), + ('S', 'skin'), ) individual = models.ForeignKey('Individual', on_delete=models.PROTECT) diff --git a/seqr/utils/search/utils.py b/seqr/utils/search/utils.py index 7bef21e164..2b52f3fe96 100644 --- a/seqr/utils/search/utils.py +++ b/seqr/utils/search/utils.py @@ -11,7 +11,7 @@ from seqr.utils.logging_utils import SeqrLogger from seqr.utils.redis_utils import safe_redis_get_json, safe_redis_get_wildcard_json, safe_redis_set_json from clickhouse_search.constants import XPOS_SORT_KEY, PRIORITIZED_GENE_SORT, RECESSIVE, COMPOUND_HET, \ - SV_ANNOTATION_TYPES, X_LINKED_RECESSIVE, PATHOGENICTY_SORT_KEY, PATHOGENICTY_HGMD_SORT_KEY + PATHOGENICTY_SORT_KEY, PATHOGENICTY_HGMD_SORT_KEY from seqr.utils.gene_utils import parse_locus_list_items from seqr.utils.xpos_utils import get_xpos, format_chrom from seqr.views.utils.permissions_utils import user_is_analyst @@ -46,11 +46,10 @@ def _get_search_genome_version(search_model): def get_single_variant(family, variant_id, user): parsed_variant_id = parse_variant_id(variant_id) - dataset_type = _variant_id_dataset_type(parsed_variant_id) parsed_variant_ids = [parsed_variant_id] if parsed_variant_id else None search = {'parsed_variant_ids': parsed_variant_ids, 'variant_ids': [variant_id]} genome_version = family.project.genome_version - variants = get_clickhouse_variants([family], [dataset_type], search, user, genome_version) + variants = get_clickhouse_variants([family], search, user, genome_version) if not variants: raise InvalidSearchException('Variant {} not found'.format(variant_id)) return format_clickhouse_results(variants, genome_version)[0] @@ -68,22 +67,12 @@ def variant_lookup(user, variant_id, genome_version, sample_type=None, affected_ return variants parsed_variant_id = parse_variant_id(variant_id) - dataset_type = _variant_id_dataset_type(parsed_variant_id) - _validate_dataset_type_genome_version(dataset_type, sample_type, genome_version) - - variants = clickhouse_variant_lookup(user, variant_id, parsed_variant_id, dataset_type, sample_type, genome_version, affected_only, hom_only) + variants = clickhouse_variant_lookup(user, variant_id, parsed_variant_id, sample_type, genome_version, affected_only, hom_only) safe_redis_set_json(cache_key, variants, expire=timedelta(weeks=2)) return variants -def _validate_dataset_type_genome_version(dataset_type, sample_type, genome_version): - if genome_version == GENOME_VERSION_GRCh37 and dataset_type != Sample.DATASET_TYPE_VARIANT_CALLS: - raise InvalidSearchException(f'{dataset_type} variants are not available for GRCh37') - if dataset_type == Sample.DATASET_TYPE_SV_CALLS and not sample_type: - raise InvalidSearchException('Sample type must be specified to look up a structural variant') - - def _get_search_cache_key(search_model, sort=None): return 'search_results__{}__{}'.format(search_model.guid, sort or XPOS_SORT_KEY) @@ -106,7 +95,7 @@ def _get_previous_search_results(search_model, sort): previous_search_results = safe_redis_get_wildcard_json(wildcard_cache_key) if previous_search_results and sort: previous_search_results = get_sorted_search_results( - previous_search_results, sort, family_guid=search_model.families.first().guid, + previous_search_results, sort, families=search_model.families.all(), ) cache_key = _get_search_cache_key(search_model, sort=sort) safe_redis_set_json(cache_key, previous_search_results, expire=timedelta(weeks=2)) @@ -134,10 +123,9 @@ def _query_variants(search_model, user, sort=None): _validate_sort(sort, families) parsed_search = _parse_search(search, genome_version, user) - dataset_types, secondary_dataset_types = _search_dataset_type(parsed_search, genome_version) _validate_search(parsed_search, families) - results = get_clickhouse_variants(families, dataset_types, parsed_search, user, genome_version, sort=sort, secondary_dataset_types=secondary_dataset_types) + results = get_clickhouse_variants(families, parsed_search, user, genome_version, sort=sort) cache_key = _get_search_cache_key(search_model, sort=sort) safe_redis_set_json(cache_key, results, expire=timedelta(weeks=2)) @@ -268,64 +256,6 @@ def _validate_sort(sort, families): raise InvalidSearchException('Phenotype sort is only supported for single-family search.') -def _search_dataset_type(search, genome_version): - parsed_variant_ids = search.get('parsed_variant_ids') - secondary_dataset_types = None - if parsed_variant_ids: - dataset_types = _chromosome_filter_dataset_types([vid[0] for vid in parsed_variant_ids]) - else: - chroms = [gene[f'chromGrch{genome_version}'] for gene in (search.get('genes') or {}).values()] + [ - interval['chrom'] for interval in (search.get('intervals') or []) - ] if not search.get('exclude_locations') else None - dataset_types = _annotation_dataset_type(search.get('annotations'), chroms, pathogenicity=search.get('pathogenicity'), exclude_svs=search.pop('exclude_svs', False)) - secondary_dataset_types = _annotation_dataset_type(search['annotations_secondary'], chroms) if search.get('annotations_secondary') else None - if secondary_dataset_types and len(dataset_types or []) == 1 and dataset_types[0] == Sample.DATASET_TYPE_SV_CALLS: - secondary_dataset_types = [dt for dt in secondary_dataset_types if dt != Sample.DATASET_TYPE_MITO_CALLS] - - if search.get('inheritance_mode') == X_LINKED_RECESSIVE: - if not dataset_types: - dataset_types = [Sample.DATASET_TYPE_VARIANT_CALLS, Sample.DATASET_TYPE_SV_CALLS] - elif Sample.DATASET_TYPE_MITO_CALLS in dataset_types: - dataset_types.remove(Sample.DATASET_TYPE_MITO_CALLS) - - return dataset_types, secondary_dataset_types - - -def _variant_id_dataset_type(parsed_variant_id): - if not parsed_variant_id: - return Sample.DATASET_TYPE_SV_CALLS - if parsed_variant_id[0].replace('chr', '').startswith('M'): - return Sample.DATASET_TYPE_MITO_CALLS - return Sample.DATASET_TYPE_VARIANT_CALLS - - -def _chromosome_filter_dataset_types(chroms): - has_mito = [chrom for chrom in chroms if chrom.replace('chr', '').startswith('M')] - if len(has_mito) == len(chroms): - return [Sample.DATASET_TYPE_MITO_CALLS] - dataset_types = [Sample.DATASET_TYPE_VARIANT_CALLS] - if has_mito: - dataset_types.append(Sample.DATASET_TYPE_MITO_CALLS) - return dataset_types - - -def _annotation_dataset_type(annotations, chroms, pathogenicity=None, exclude_svs=False): - if not (annotations or chroms): - return [Sample.DATASET_TYPE_VARIANT_CALLS, Sample.DATASET_TYPE_MITO_CALLS] if (pathogenicity or exclude_svs) else None - - annotation_types = set((annotations or {}).keys()) - if annotations and annotation_types.issubset(SV_ANNOTATION_TYPES) and not pathogenicity: - return [Sample.DATASET_TYPE_SV_CALLS] - - dataset_types = _chromosome_filter_dataset_types(chroms) if chroms else [Sample.DATASET_TYPE_VARIANT_CALLS, Sample.DATASET_TYPE_MITO_CALLS] - - no_svs = exclude_svs or (annotations and annotation_types.isdisjoint(SV_ANNOTATION_TYPES)) - if not no_svs: - dataset_types.append(Sample.DATASET_TYPE_SV_CALLS) - - return dataset_types - - def _parse_inheritance(search): inheritance = search.pop('inheritance') inheritance_mode = inheritance.get('mode') diff --git a/seqr/views/apis/data_manager_api.py b/seqr/views/apis/data_manager_api.py index 821bff7af3..6f0fc97ae9 100644 --- a/seqr/views/apis/data_manager_api.py +++ b/seqr/views/apis/data_manager_api.py @@ -37,7 +37,9 @@ 'muscle': 'Muscle', 'airway_cultured_epithelium': 'Nasal Epithelium', 'brain': 'Brain', - 'lymphocytes': 'LCL Pellet' + 'lymphocytes': 'LCL Pellet', + 'liver': 'Liver', + 'skin': 'Skin Epidermis', } TISSUE_TYPE_MAP = { AIRTABLE_TISSUE_TYPE_MAP[name]: type diff --git a/ui/package-lock.json b/ui/package-lock.json index 022be7b520..c1cab66e44 100644 --- a/ui/package-lock.json +++ b/ui/package-lock.json @@ -13093,9 +13093,9 @@ } }, "node_modules/node-forge": { - "version": "1.3.2", - "resolved": "https://registry.npmjs.org/node-forge/-/node-forge-1.3.2.tgz", - "integrity": "sha512-6xKiQ+cph9KImrRh0VsjH2d8/GXA4FIMlgU4B757iI1ApvcyA9VlouP0yZJha01V+huImO+kKMU7ih+2+E14fw==", + "version": "1.4.0", + "resolved": "https://registry.npmjs.org/node-forge/-/node-forge-1.4.0.tgz", + "integrity": "sha512-LarFH0+6VfriEhqMMcLX2F7SwSXeWwnEAJEsYm5QKWchiVYVvJyV9v7UDvUv+w5HO23ZpQTXDv/GxdDdMyOuoQ==", "dev": true, "engines": { "node": ">= 6.13.0" @@ -29518,9 +29518,9 @@ } }, "node-forge": { - "version": "1.3.2", - "resolved": "https://registry.npmjs.org/node-forge/-/node-forge-1.3.2.tgz", - "integrity": "sha512-6xKiQ+cph9KImrRh0VsjH2d8/GXA4FIMlgU4B757iI1ApvcyA9VlouP0yZJha01V+huImO+kKMU7ih+2+E14fw==", + "version": "1.4.0", + "resolved": "https://registry.npmjs.org/node-forge/-/node-forge-1.4.0.tgz", + "integrity": "sha512-LarFH0+6VfriEhqMMcLX2F7SwSXeWwnEAJEsYm5QKWchiVYVvJyV9v7UDvUv+w5HO23ZpQTXDv/GxdDdMyOuoQ==", "dev": true }, "node-int64": { diff --git a/ui/pages/Project/selectors.js b/ui/pages/Project/selectors.js index 172b8ad08b..7fb9dafb6b 100644 --- a/ui/pages/Project/selectors.js +++ b/ui/pages/Project/selectors.js @@ -136,7 +136,7 @@ export const getProjectAnalysisGroupFamilySizeHistogram = createSelector( export const getProjectAnalysisGroupDataLoadedFamilySizeHistogram = createSelector( getProjectAnalysisGroupFamiliesByGuid, getActiveDatasetsByFamily, - (familiesByGuid, datasetsByFamily) => getFamilySizeHistogram(Object.values(datasetsByFamily).map(((family) => { + (familiesByGuid, datasetsByFamily) => getFamilySizeHistogram(Object.values(familiesByGuid).map(((family) => { const sampleIndividuals = new Set((datasetsByFamily[family.familyGuid] || []).map(sample => sample.individualGuid)) const hasSampleParents = (family.parents || []).reduce( (acc, { individualGuid, maternalGuid, paternalGuid }) => { diff --git a/ui/pages/Search/VariantLookup.jsx b/ui/pages/Search/VariantLookup.jsx index ccbb4f6982..e47264c8ba 100644 --- a/ui/pages/Search/VariantLookup.jsx +++ b/ui/pages/Search/VariantLookup.jsx @@ -73,7 +73,7 @@ const LOOKUP_FIELDS = [ ] const validateAnnotations = (value, { search }) => ( - value || Object.values(search.annotations || {}).some(val => val.length) ? undefined : 'At least one consequence filter is required' + value || Object.values(search?.annotations || {}).some(val => val.length) ? undefined : 'At least one consequence filter is required' ) const CONSEQUENCE_FIELDS = [