Enforce a 1 to 1 relationship between incoming score rows and outgoing annotated mappings

Prior to this change, it was possible for some score rows to generate valid mappings with other score rows not creating a mapped variant. This had some negative downstream consequences, which will be remedied by ensuring that if any variant receives a mapped variant, all variants receive a mapped variant.

Author: bencap

src/api/routers/map.py

@@ -68,6 +68,7 @@ async def map_scoreset(urn: str, store_path: Path | None = None) -> JSONResponse
                 error_message="Score set contains no variants to map",
             ).model_dump(exclude_none=True)
         )
+    total_score_records = sum(len(v) for v in records.values())
 
     try:
         alignment_results = build_alignment_result(metadata, True)
@@ -120,15 +121,27 @@ async def map_scoreset(urn: str, store_path: Path | None = None) -> JSONResponse
                 metadata=metadata, error_message=str(e).strip("'")
             ).model_dump(exclude_none=True)
         )
-    if not vrs_results or all(
-        mapping_result is None for mapping_result in vrs_results.values()
-    ):
+
+    nonetype_vrs_results = [
+        result is None
+        for target_gene in vrs_results
+        for result in vrs_results[target_gene]
+    ]
+
+    if not vrs_results or all(nonetype_vrs_results):
         return JSONResponse(
             content=ScoresetMapping(
                 metadata=metadata,
                 error_message="No variant mappings available for this score set",
             ).model_dump(exclude_none=True)
         )
+    if any(nonetype_vrs_results):
+        return JSONResponse(
+            content=ScoresetMapping(
+                metadata=metadata,
+                error_message="Some variants generated vrs results, but not all. If any variants were mapped, all should have been.",
+            ).model_dump(exclude_none=True)
+        )
 
     annotated_vrs_results = {}
     try:
@@ -146,15 +159,27 @@ async def map_scoreset(urn: str, store_path: Path | None = None) -> JSONResponse
                 metadata=metadata, error_message=str(e).strip("'")
             ).model_dump(exclude_none=True)
         )
-    if not annotated_vrs_results or all(
-        mapping_result is None for mapping_result in annotated_vrs_results.values()
-    ):
+
+    nonetype_annotated_vrs_results = [
+        result is None
+        for target_gene in annotated_vrs_results
+        for result in annotated_vrs_results[target_gene]
+    ]
+
+    if not annotated_vrs_results or all(nonetype_annotated_vrs_results):
         return JSONResponse(
             content=ScoresetMapping(
                 metadata=metadata,
                 error_message="No annotated variant mappings available for this score set",
             ).model_dump(exclude_none=True)
         )
+    if any(nonetype_annotated_vrs_results):
+        return JSONResponse(
+            content=ScoresetMapping(
+                metadata=metadata,
+                error_message="Some variants generated annotated vrs results, but not all. If any variants were annotated, all should have been.",
+            ).model_dump(exclude_none=True)
+        )
 
     try:
         raw_metadata = get_raw_scoreset_metadata(urn, store_path)
@@ -235,6 +260,14 @@ async def map_scoreset(urn: str, store_path: Path | None = None) -> JSONResponse
             ).model_dump(exclude_none=True)
         )
 
+    if len(mapped_scores) != total_score_records:
+        return JSONResponse(
+            content=ScoresetMapping(
+                metadata=metadata,
+                error_message=f"Mismatch between number of mapped scores ({len(mapped_scores)}) and total score records ({total_score_records}). This is unexpected and indicates an issue with the mapping process.",
+            ).model_dump(exclude_none=True)
+        )
+
     return JSONResponse(
         content=ScoresetMapping(
             metadata=raw_metadata,

src/dcd_mapping/annotate.py

@@ -266,13 +266,16 @@ def _annotate_allele_mapping(
             sequence_id = f"ga4gh:{mapped_score.post_mapped.location.sequenceReference.refgetAccession}"
             accession = get_chromosome_identifier_from_vrs_id(sequence_id)
             if accession is None:
-                raise ValueError
-            if accession.startswith("refseq:"):
+                accession = None
+                mapped_score.error_message = "Could not determine accession for this annotation. No allele expression is available."
+            elif accession.startswith("refseq:"):
                 accession = accession[7:]
         else:
             if tx_results is None or isinstance(tx_results, TxSelectError):
-                raise ValueError  # impossible by definition
-            accession = tx_results.np
+                accession = None
+                mapped_score.error_message = "Could not determine accession for this annotation. No allele expression is available."
+            else:
+                accession = tx_results.np
 
         sr = get_seqrepo()
         loc = mapped_score.post_mapped.location
@@ -281,8 +284,9 @@ def _annotate_allele_mapping(
         post_mapped.extensions = [
             Extension(type="Extension", name="vrs_ref_allele_seq", value=ref)
         ]
-        hgvs_string, syntax = _get_hgvs_string(post_mapped, accession)
-        post_mapped.expressions = [Expression(syntax=syntax, value=hgvs_string)]
+        if accession:
+            hgvs_string, syntax = _get_hgvs_string(post_mapped, accession)
+            post_mapped.expressions = [Expression(syntax=syntax, value=hgvs_string)]
 
     if vrs_version == VrsVersion.V_1_3:
         pre_mapped = _allele_to_vod(pre_mapped)
@@ -295,9 +299,7 @@ def _annotate_allele_mapping(
         mavedb_id=mapped_score.accession_id,
         score=float(mapped_score.score) if mapped_score.score else None,
         annotation_layer=mapped_score.annotation_layer,
-        error_message=mapped_score.error_message
-        if mapped_score.error_message
-        else None,  # TODO might not need if statement here
+        error_message=mapped_score.error_message,
     )
 
 
@@ -321,13 +323,17 @@ def _annotate_haplotype_mapping(
             sequence_id = f"ga4gh:{post_mapped.members[0].location.sequenceReference.refgetAccession}"
             accession = get_chromosome_identifier_from_vrs_id(sequence_id)
             if accession is None:
-                raise ValueError
-            if accession.startswith("refseq:"):
+                accession = None
+                mapped_score.error_message = "Could not determine accession for this annotation. No allele expression is available."
+            elif accession.startswith("refseq:"):
                 accession = accession[7:]
         else:
             if tx_results is None or isinstance(tx_results, TxSelectError):
-                raise ValueError  # impossible by definition
-            accession = tx_results.np
+                # impossible by definition
+                accession = None
+                mapped_score.error_message = "Could not determine accession for this annotation. No allele expression is available."
+            else:
+                accession = tx_results.np
 
         sr = get_seqrepo()
         for allele in post_mapped.members:
@@ -337,8 +343,9 @@ def _annotate_haplotype_mapping(
             allele.extensions = [
                 Extension(type="Extension", name="vrs_ref_allele_seq", value=ref)
             ]
-            hgvs, syntax = _get_hgvs_string(allele, accession)
-            allele.expressions = [Expression(syntax=syntax, value=hgvs)]
+            if accession:
+                hgvs, syntax = _get_hgvs_string(allele, accession)
+                allele.expressions = [Expression(syntax=syntax, value=hgvs)]
 
     if vrs_version == VrsVersion.V_1_3:
         pre_mapped = _haplotype_to_haplotype_1_3(pre_mapped)
@@ -351,9 +358,7 @@ def _annotate_haplotype_mapping(
         mavedb_id=mapped_score.accession_id,
         score=float(mapped_score.score) if mapped_score.score is not None else None,
         annotation_layer=mapped_score.annotation_layer,
-        error_message=mapped_score.error_message
-        if mapped_score.error_message
-        else None,  # TODO might not need if statement here
+        error_message=mapped_score.error_message,
     )
 
 
@@ -388,6 +393,7 @@ def annotate(
                 ScoreAnnotationWithLayer(
                     mavedb_id=mapped_score.accession_id,
                     score=float(mapped_score.score) if mapped_score.score else None,
+                    vrs_version=vrs_version,
                     error_message=mapped_score.error_message,
                 )
             )
@@ -410,8 +416,16 @@ def annotate(
                 )
             )
         else:
-            # TODO how to combine this error message with other potential error messages?
-            ValueError("inconsistent variant structure")
+            score_annotations.append(
+                ScoreAnnotationWithLayer(
+                    pre_mapped=mapped_score.pre_mapped,
+                    post_mapped=mapped_score.post_mapped,
+                    vrs_version=vrs_version,
+                    mavedb_id=mapped_score.accession_id,
+                    score=float(mapped_score.score) if mapped_score.score else None,
+                    error_message=f"Multiple issues with annotation: Inconsistent variant structure (Allele and Haplotype mix).{' ' + mapped_score.error_message if mapped_score.error_message else ''}",
+                )
+            )
 
     return score_annotations
 
@@ -519,11 +533,6 @@ def _set_scoreset_layer(
     expressions. If genomic expressions are available, that's what we'd like to use.
     This function tells us how to filter the final annotation objects.
     """
-    if urn.startswith("urn:mavedb:00000097"):
-        _logger.debug(
-            "Manually selecting protein annotation for scores from urn:mavedb:00000097"
-        )
-        return AnnotationLayer.PROTEIN
     for mapping in mappings:
         if mapping.annotation_layer == AnnotationLayer.GENOMIC:
             return AnnotationLayer.GENOMIC

src/dcd_mapping/transcripts.py

@@ -54,7 +54,10 @@ async def _get_compatible_transcripts(
     chromosome = get_chromosome_identifier(aligned_chrom)
     gene_symbol = get_gene_symbol(target_gene)
     if not gene_symbol:
-        raise TxSelectError
+        msg = (
+            f"Unable to find gene symbol for target gene {target_gene.target_gene_name}"
+        )
+        raise TxSelectError(msg)
     transcript_matches = []
     for hit_range in align_result.hit_subranges:
         matches_list = await get_transcripts(
@@ -179,7 +182,8 @@ async def _select_protein_reference(
         if not best_tx:
             best_tx = await _get_longest_compatible_transcript(common_transcripts)
         if not best_tx:
-            raise TxSelectError
+            msg = f"Unable to find matching MANE transcripts for target gene {target_gene.target_gene_name}"
+            raise TxSelectError(msg)
         ref_sequence = get_sequence(best_tx.refseq_prot)
         nm_accession = best_tx.refseq_nuc
         np_accession = best_tx.refseq_prot
@@ -323,19 +327,6 @@ async def select_transcript(
     :param align_result: alignment results
     :return: Transcript description (accession ID, offset, selected sequence, etc)
     """
-    if scoreset_urn.startswith("urn:mavedb:00000097"):
-        _logger.info(
-            "Score sets in urn:mavedb:00000097 are already expressed in full HGVS strings -- using predefined results because additional hard-coding is unnecessary"
-        )
-        return TxSelectResult(
-            nm="NM_007294.3",
-            np="NP_009225.1",
-            start=0,
-            is_full_match=False,
-            transcript_mode=TranscriptPriority.MANE_SELECT,
-            sequence=_get_protein_sequence(target_gene.target_sequence),
-        )
-
     if target_gene.target_gene_category != TargetType.PROTEIN_CODING:
         _logger.debug(
             "%s is regulatory/noncoding -- skipping transcript selection",
@@ -366,7 +357,9 @@ async def select_transcripts(
         * Dict where keys are target labels and values are alignment result objects
     :return: dict where keys are target labels and values are objects describing selected transcript (accession ID, offset, selected sequence, etc)
     """
-    selected_transcripts = {}
+    selected_transcripts: dict[
+        str, TxSelectResult | TxSelectError | KeyError | None
+    ] = {}
     for target_gene in scoreset_metadata.target_genes:
         if scoreset_metadata.target_genes[target_gene].target_accession_id:
             # for accession-based targets, create tx select objects for protein sequence accessions only

src/dcd_mapping/vrs_map.py

@@ -453,8 +453,11 @@ def _map_genomic(
             # if the sequence id starts with SQ, it is a target sequence which is in the ga4gh namespace
             namespace = "ga4gh"
         else:
-            msg = f"Namespace could not be inferred from sequence: {sequence_id}"
-            raise ValueError(msg)
+            return MappedScore(
+                accession_id=row.accession,
+                score=row.score,
+                error_message=f"Namespace could not be inferred from sequence: {sequence_id}",
+            )
 
     if (
         row.hgvs_nt in {"_wt", "_sy", "="}
@@ -609,8 +612,11 @@ def _map_genomic(
                 error_message=str(e),
             )
     else:
-        msg = f"Reference sequence namespace not supported: {namespace}"
-        raise ValueError(msg)
+        return MappedScore(
+            accession_id=row.accession,
+            score=row.score,
+            error_message=f"Reference sequence namespace not supported: {namespace}",
+        )
 
     return MappedScore(
         accession_id=row.accession,
@@ -783,7 +789,14 @@ def _map_accession(
     variations: list[MappedScore] = []
     sequence_id = metadata.target_accession_id
     if sequence_id is None:
-        raise ValueError
+        return [
+            MappedScore(
+                accession_id=row.accession,
+                score=row.score,
+                error_message="Could not generate mapped allele objects. No sequence id was provided.",
+            )
+            for row in records
+        ]
 
     store_accession(sequence_id)
 
@@ -802,8 +815,14 @@ def _map_accession(
             hgvs_nt_mappings = _map_genomic(row, sequence_id, align_result)
             variations.append(hgvs_nt_mappings)
     else:
-        msg = f"Unrecognized accession prefix for accession id {metadata.target_accession_id}"
-        raise ValueError(msg)
+        [
+            MappedScore(
+                accession_id=row.accession,
+                score=row.score,
+                error_message=f"Unrecognized accession prefix for accession id {metadata.target_accession_id}",
+            )
+            for row in records
+        ]
 
     return variations
 
@@ -887,7 +906,7 @@ def vrs_map(
     records: list[ScoreRow],
     transcript: TxSelectResult | TxSelectError | None = None,
     silent: bool = True,
-) -> list[MappedScore] | None:
+) -> list[MappedScore]:
     """Given a description of a MAVE scoreset and an aligned transcript, generate
     the corresponding VRS objects.