feat(grz-pydantic-models): Add a thresholds model and implement check for oncomine panel submissions (#469)

Fixes #462

Author: Hoeze

packages/grz-common/src/grz_common/workers/submission.py

@@ -11,18 +11,16 @@
 from os import PathLike
 from pathlib import Path
 
+import grz_pydantic_models.submission.thresholds as thresholds_model
 from grz_pydantic_models.submission.metadata import get_accepted_versions
 from grz_pydantic_models.submission.metadata.v1 import (
     ChecksumType,
     File,
     FileType,
     GrzSubmissionMetadata,
-    LibraryType,
     ReadOrder,
     SequenceData,
-    SequenceSubtype,
     SequencingLayout,
-    load_thresholds,
 )
 from grz_pydantic_models.submission.metadata.v1 import File as SubmissionFileMetadata
 from pydantic import ValidationError
@@ -236,15 +234,8 @@ def should_check_file(file_path: Path, file_metadata: File) -> bool:
                 if not lab_data.sequence_data:
                     continue
 
-                threshold_definitions = load_thresholds()
-                thresholds = threshold_definitions[
-                    (
-                        self.metadata.content.submission.genomic_study_subtype,
-                        lab_data.library_type,
-                        lab_data.sequence_subtype,
-                    )
-                ]
-                mean_read_length_threshold = thresholds.get("meanReadLength", 0)
+                thresholds = self.metadata.content.determine_thresholds_for(donor, lab_data)
+                mean_read_length_threshold = thresholds.mean_read_length
 
                 sequence_data = lab_data.sequence_data
                 fastq_files = [f for f in sequence_data.files if f.file_type == FileType.fastq]
@@ -501,22 +492,26 @@ def find_bam_files(sequence_data: SequenceData) -> list[File]:
             return [f for f in sequence_data.files if f.file_type == FileType.bam]
 
         for donor in self.metadata.content.donors:
-            for lab_data in donor.lab_data:
-                sequencing_layout = lab_data.sequencing_layout
-                sequence_data = lab_data.sequence_data
+            for lab_datum in donor.lab_data:
+                sequencing_layout = lab_datum.sequencing_layout
+                sequence_data = lab_datum.sequence_data
                 # find all FASTQ files
                 fastq_files = find_fastq_files(sequence_data) if sequence_data else []
                 bam_files = find_bam_files(sequence_data) if sequence_data else []
 
-                if not lab_data.library_type.endswith("_lr"):
+                if not lab_datum.library_type.endswith("_lr"):
                     match sequencing_layout:
                         case SequencingLayout.single_end | SequencingLayout.reverse | SequencingLayout.other:
                             yield from self._validate_single_end_fallback(
-                                fastq_files, progress_logger, lab_data.library_type, lab_data.sequence_subtype
+                                fastq_files,
+                                progress_logger,
+                                self.metadata.content.determine_thresholds_for(donor, lab_datum),
                             )
                         case SequencingLayout.paired_end:
                             yield from self._validate_paired_end_fallback(
-                                fastq_files, progress_logger, lab_data.library_type, lab_data.sequence_subtype
+                                fastq_files,
+                                progress_logger,
+                                self.metadata.content.determine_thresholds_for(donor, lab_datum),
                             )
                 yield from self._validate_bams_fallback(bam_files, progress_logger)
 
@@ -559,22 +554,13 @@ def _validate_single_end_fallback(
         self,
         fastq_files: list[File],
         progress_logger: FileProgressLogger[ValidationState],
-        library_type: LibraryType,
-        sequence_subtype: SequenceSubtype,
+        thresholds: thresholds_model.Thresholds,
     ) -> Generator[str, None, None]:
         def validate_file(local_file_path, file_metadata: SubmissionFileMetadata) -> ValidationState:
             self.__log.debug("Validating '%s'...", str(local_file_path))
 
             # validate the file
-            threshold_definitions = load_thresholds()
-            thresholds = threshold_definitions[
-                (
-                    self.metadata.content.submission.genomic_study_subtype,
-                    library_type,
-                    sequence_subtype,
-                )
-            ]
-            mean_read_length_threshold = thresholds["meanReadLength"]
+            mean_read_length_threshold = thresholds.mean_read_length
             errors = list(
                 validate_single_end_reads(local_file_path, mean_read_length_threshold=mean_read_length_threshold)
             )
@@ -596,18 +582,9 @@ def _validate_paired_end_fallback(
         self,
         fastq_files: list[File],
         progress_logger: FileProgressLogger[ValidationState],
-        library_type: LibraryType,
-        sequence_subtype: SequenceSubtype,
+        thresholds: thresholds_model.Thresholds,
     ) -> Generator[str, None, None]:
-        threshold_definitions = load_thresholds()
-        thresholds = threshold_definitions[
-            (
-                self.metadata.content.submission.genomic_study_subtype,
-                library_type,
-                sequence_subtype,
-            )
-        ]
-        mean_read_length_threshold = thresholds["meanReadLength"]
+        mean_read_length_threshold = thresholds.mean_read_length
         key = lambda f: (f.flowcell_id, f.lane_id)
         fastq_files.sort(key=key)
         for _key, group in groupby(fastq_files, key):

packages/grz-pydantic-models/src/grz_pydantic_models/submission/metadata/v1.py

@@ -6,11 +6,12 @@
 import re
 from datetime import date
 from enum import StrEnum
+from functools import cache
 from importlib.resources import files
 from itertools import groupby
 from operator import attrgetter
 from pathlib import Path, PurePosixPath
-from typing import Annotated, Any, Self
+from typing import Annotated, Self
 
 from pydantic import (
     AfterValidator,
@@ -25,6 +26,7 @@
 
 from ...common import StrictBaseModel
 from ...mii.consent import Consent, ProvisionType
+from .. import thresholds as thresholds_model
 from .versioning import Version
 
 SCHEMA_URL_PATTERN = r"https://raw\.githubusercontent\.com/BfArM-MVH/MVGenomseq/refs/tags/v([0-9]+)\.([0-9]+)(?:\.([0-9]+))?/GRZ/grz-schema\.json"
@@ -1290,40 +1292,81 @@ def check_duplicate_file_paths(self):
 
         return self
 
+    def determine_thresholds_for(
+        self,
+        donor: Donor,
+        lab_datum: LabDatum,
+    ) -> thresholds_model.Thresholds:
+        """
+        Determine the thresholds for a given donor and lab datum.
+
+        :param donor: The donor for which to determine the thresholds.
+        :param lab_datum: The lab datum for which to determine the thresholds.
+        :raise MissingThresholdsException: If no thresholds are found for the given combination.
+        :return: The thresholds for the given donor and lab datum.
+        """
+        threshold_definitions = _load_thresholds()
+
+        genomic_study_subtype = self.submission.genomic_study_subtype
+        is_oncomine_panel = "oncomine" in lab_datum.kit_name.lower()
+        key = (genomic_study_subtype, lab_datum.library_type, lab_datum.sequence_subtype, is_oncomine_panel)
+
+        thresholds = threshold_definitions.get(key)
+        if thresholds is None:
+            allowed_combinations_list = sorted(list(threshold_definitions.keys()))
+            allowed_combinations = "\n".join([f"  - {combination}" for combination in allowed_combinations_list])
+            names = (
+                "submission.genomicStudySubtype",
+                "labData.libraryType",
+                "labData.sequenceSubtype",
+                "isOncominePanel",
+            )
+
+            info = dict(zip(names, key, strict=True))
+
+            raise MissingThresholdsException(
+                f"No thresholds for the specified combination {info} found (donor {donor.donor_pseudonym})!\n"
+                f"Valid combinations:\n{allowed_combinations}.\n"
+                f"See https://www.bfarm.de/SharedDocs/Downloads/DE/Forschung/modellvorhaben-genomsequenzierung/Qs-durch-GRZ.pdf?__blob=publicationFile for more details.\n"
+                f"Skipping threshold validation."
+            )
+
+        return thresholds
+
+    def get_thresholds(self) -> dict[tuple[str, str], thresholds_model.Thresholds | None]:
+        """
+        Get the thresholds for all lab data in the submission.
+
+        :return: A dictionary mapping (donor pseudonym, lab data name) to the thresholds,
+                 or None if no thresholds are defined for this lab datum.
+        """
+        thresholds_dict = {}
+        for donor in self.donors:
+            for lab_datum in donor.lab_data:
+                try:
+                    thresholds = self.determine_thresholds_for(donor, lab_datum)
+                except MissingThresholdsException:
+                    # assign None if no thresholds are found
+                    thresholds = None
+
+                thresholds_dict[(donor.donor_pseudonym, lab_datum.lab_data_name)] = thresholds
+
+        return thresholds_dict
+
     @model_validator(mode="after")
     def validate_thresholds(self):
         """
         Check if the submission meets the minimum mean coverage requirements.
         """
-        threshold_definitions = load_thresholds()
-
         for donor in self.donors:
             for lab_datum in donor.lab_data:
-                key = (
-                    self.submission.genomic_study_subtype,
-                    lab_datum.library_type,
-                    lab_datum.sequence_subtype,
-                )
-                thresholds = threshold_definitions.get(key)
-                if thresholds is None:
-                    allowed_combinations = sorted(list(threshold_definitions.keys()))
-                    allowed_combinations = "\n".join([f"  - {combination}" for combination in allowed_combinations])
-                    names = (
-                        "submission.genomicStudySubtype",
-                        "labData.libraryType",
-                        "labData.sequenceSubtype",
-                    )
-                    info = dict(zip(names, key, strict=True))
-                    log.warning(
-                        f"No thresholds for the specified combination {info} found (donor {donor.donor_pseudonym})!\n"
-                        f"Valid combinations:\n{allowed_combinations}.\n"
-                        f"See https://www.bfarm.de/SharedDocs/Downloads/DE/Forschung/modellvorhaben-genomsequenzierung/Qs-durch-GRZ.pdf?__blob=publicationFile for more details.\n"
-                        f"Skipping threshold validation."
-                    )
+                try:
+                    thresholds = self.determine_thresholds_for(donor, lab_datum)
+                except MissingThresholdsException as e:
+                    log.warning(str(e))
                     continue
 
                 _check_thresholds(donor, lab_datum, thresholds)
-
         return self
 
     @model_validator(mode="after")
@@ -1345,7 +1388,7 @@ def validate_reference_genome_compatibility(self):
         return self
 
 
-def _check_thresholds(donor: Donor, lab_datum: LabDatum, thresholds: dict[str, Any]):  # noqa: C901
+def _check_thresholds(donor: Donor, lab_datum: LabDatum, thresholds: thresholds_model.Thresholds):  # noqa: C901
     def raise_if_index(message: str) -> None:
         if donor.relation == Relation.index_:
             raise ValueError(message)
@@ -1359,36 +1402,36 @@ def raise_if_index(message: str) -> None:
     lab_data_name = lab_datum.lab_data_name
     sequence_data = lab_datum.sequence_data
 
-    mean_depth_of_coverage_t = thresholds.get("meanDepthOfCoverage")
+    mean_depth_of_coverage_t = thresholds.mean_depth_of_coverage
     mean_depth_of_coverage_v = sequence_data.mean_depth_of_coverage
     if mean_depth_of_coverage_t and mean_depth_of_coverage_v < mean_depth_of_coverage_t:
         raise_if_index(
             f"Mean depth of coverage for donor '{pseudonym}', lab datum '{lab_data_name}' "
             f"below threshold: {mean_depth_of_coverage_v} < {mean_depth_of_coverage_t}"
         )
 
-    if percent_bases_above_quality_threshold_t := thresholds.get("percentBasesAboveQualityThreshold"):
-        minimum_quality_t = percent_bases_above_quality_threshold_t.get("qualityThreshold")
+    if percent_bases_above_quality_threshold_t := thresholds.percent_bases_above_quality_threshold:
+        minimum_quality_t = percent_bases_above_quality_threshold_t.quality_threshold
         minimum_quality_v = sequence_data.percent_bases_above_quality_threshold.minimum_quality
         if minimum_quality_t and (minimum_quality_t != minimum_quality_v):
             # TODO also print out genomic study subtype because that determines thresholds, but is defined at submission level
             # this should error regardless of relation because it should match the official threshold
             raise ValueError(
-                f"Expected minimumQuality '{minimum_quality_t}' for library type '{lab_datum.library_type}'"
+                f"Expected minimumQuality '{minimum_quality_t}' for library type '{lab_datum.library_type}' "
                 f"and sequence subtype '{lab_datum.sequence_subtype}'. Got '{minimum_quality_v}' instead."
             )
-        percent_t = percent_bases_above_quality_threshold_t.get("percentBasesAbove")
+        percent_t = percent_bases_above_quality_threshold_t.percent_bases_above
         percent_v = sequence_data.percent_bases_above_quality_threshold.percent
         if percent_t and (percent_v < percent_t):
             raise_if_index(
                 f"Percentage of bases above quality threshold are below threshold: {percent_v} < {percent_t}"
             )
 
-    if t := thresholds.get("targetedRegionsAboveMinCoverage"):
-        min_coverage_t = t.get("minCoverage")
+    if t := thresholds.targeted_regions_above_min_coverage:
+        min_coverage_t = t.min_coverage
         min_coverage_v = sequence_data.min_coverage
 
-        fraction_above_t = t.get("fractionAbove")
+        fraction_above_t = t.fraction_above
         fraction_above_v = sequence_data.targeted_regions_above_min_coverage
 
         if min_coverage_t and min_coverage_v != min_coverage_t:
@@ -1405,15 +1448,59 @@ def raise_if_index(message: str) -> None:
             )
 
 
-type Thresholds = dict[tuple[str, str, str], dict[str, Any]]
+class MissingThresholdsException(ValueError):
+    """Exception raised when no thresholds are found for a given combination."""
+
+    pass
+
 
+# Dictionary type for thresholds lookup:
+# Keys are (GenomicStudySubtype, LibraryType, SequenceSubtype, is_oncomine_panel)
+_ThresholdsDict = dict[
+    tuple[GenomicStudySubtype, LibraryType, SequenceSubtype, bool],
+    thresholds_model.Thresholds,
+]
 
-def load_thresholds() -> Thresholds:
-    threshold_definitions = json.load(
+
+class _ThresholdEntry(StrictBaseModel):
+    """Threshold configuration for a particular sequencing setup."""
+
+    model_config = ConfigDict(populate_by_name=True)
+
+    library_type: LibraryType = Field(alias="libraryType")
+    sequence_subtype: SequenceSubtype = Field(alias="sequenceSubtype")
+    genomic_study_subtype: GenomicStudySubtype = Field(alias="genomicStudySubtype")
+    thresholds: thresholds_model.Thresholds
+
+
+@cache
+def _load_thresholds() -> _ThresholdsDict:
+    threshold_definitions_json = json.load(
         files("grz_pydantic_models").joinpath("resources", "thresholds.json").open("r", encoding="utf-8")
     )
-    threshold_definitions = {
-        (d["genomicStudySubtype"], d["libraryType"], d["sequenceSubtype"]): d["thresholds"]
+    try:
+        threshold_definitions: list[_ThresholdEntry] = [
+            _ThresholdEntry.model_validate(item) for item in threshold_definitions_json
+        ]
+    except ValidationError as err:
+        raise RuntimeError("Invalid threshold definitions JSON") from err
+
+    threshold_definitions_dict = {
+        (d.genomic_study_subtype, d.library_type, d.sequence_subtype, False): d.thresholds
         for d in threshold_definitions
     }
-    return threshold_definitions
+    # add thresholds for oncomine panels
+    for d in threshold_definitions:
+        if d.library_type == LibraryType.panel:
+            thresholds = d.thresholds.model_copy(deep=True)
+            # we have different thresholds_model.percentBasesAboveQualityThreshold for oncomine panels:
+            # percentage of bases ≥ Q20 should be at least 70%
+            thresholds.percent_bases_above_quality_threshold = thresholds_model.PercentBasesAboveQualityThreshold(
+                qualityThreshold=20,
+                percentBasesAbove=70,
+            )
+
+            # add new entry for oncomine panels
+            threshold_definitions_dict[(d.genomic_study_subtype, d.library_type, d.sequence_subtype, True)] = thresholds
+
+    return threshold_definitions_dict

packages/grz-pydantic-models/src/grz_pydantic_models/submission/thresholds/__init__.py

@@ -0,0 +1 @@
+from .v1 import *  # noqa: F403