Validate imputed sex results against the callset. (#780)

* scaffold

* Add validation of imputed sex file

* add args

* ruff

* Update v03_pipeline/lib/misc/validation_test.py

Co-authored-by: Daniel Marten <78616802+matren395@users.noreply.github.com>

* Update v03_pipeline/lib/misc/validation_test.py

Co-authored-by: Daniel Marten <78616802+matren395@users.noreply.github.com>

* Removing support for Unknown

* More explicit

* format

* lint

---------

Co-authored-by: Daniel Marten <78616802+matren395@users.noreply.github.com>

Author: bpblanken

v03_pipeline/lib/misc/io.py

@@ -5,7 +5,6 @@
 import hail as hl
 
 from v03_pipeline.lib.misc.gcnv import parse_gcnv_genes
-from v03_pipeline.lib.misc.validation import validate_ambiguous_sex
 from v03_pipeline.lib.model import DatasetType, Env, ReferenceGenome, Sex
 
 BIALLELIC = 2
@@ -168,12 +167,15 @@ def import_imputed_sex(imputed_sex_path: str) -> hl.Table:
             hl.case()
             .when(ht.predicted_sex == FEMALE, Sex.FEMALE.value)
             .when(ht.predicted_sex == MALE, Sex.MALE.value)
-            .default(Sex.UNKNOWN.value)
+            .or_error(
+                hl.format(
+                    'Found unexpected value %s in imputed sex file',
+                    ht.predicted_sex,
+                ),
+            )
         ),
     )
-    ht = ht.key_by(ht.s)
-    validate_ambiguous_sex(ht)
-    return ht
+    return ht.key_by(ht.s)
 
 
 def import_remap(remap_path: str) -> hl.Table:

v03_pipeline/lib/misc/io_test.py

@@ -9,6 +9,9 @@
 )
 
 TEST_IMPUTED_SEX = 'v03_pipeline/var/test/sex_check/test_imputed_sex.tsv'
+TEST_IMPUTED_SEX_UNEXPECTED_VALUE = (
+    'v03_pipeline/var/test/sex_check/test_imputed_sex_unexpected_value.tsv'
+)
 TEST_MITO_MT = 'v03_pipeline/var/test/callsets/mito_1.mt'
 TEST_SV_VCF = 'v03_pipeline/var/test/callsets/sv_1.vcf'
 
@@ -35,3 +38,11 @@ def test_import_imputed_sex(self) -> None:
                 hl.Struct(s='abc_3', predicted_sex='M'),
             ],
         )
+
+    def test_import_imputed_sex_unexpected_value(self) -> None:
+        ht = import_imputed_sex(TEST_IMPUTED_SEX_UNEXPECTED_VALUE)
+        self.assertRaisesRegex(
+            hl.utils.java.HailUserError,
+            'Found unexpected value Unknown in imputed sex file',
+            ht.collect,
+        )

v03_pipeline/lib/misc/validation.py

@@ -11,15 +11,6 @@ class SeqrValidationError(Exception):
     pass
 
 
-def validate_ambiguous_sex(ht: hl.Table) -> None:
-    ambiguous_perc = ht.aggregate(
-        hl.agg.fraction(ht.predicted_sex == Sex.UNKNOWN.value),
-    )
-    if ambiguous_perc > AMBIGUOUS_THRESHOLD_PERC:
-        msg = f'{ambiguous_perc:.2%} of samples identified as ambiguous.  Please contact the methods team to investigate the callset.'
-        raise ValueError(msg)
-
-
 def validate_no_duplicate_variants(
     mt: hl.MatrixTable,
 ) -> None:
@@ -56,6 +47,30 @@ def validate_expected_contig_frequency(
             raise SeqrValidationError(msg)
 
 
+def validate_imputed_sex_ploidy(
+    mt: hl.MatrixTable,
+    sex_check_ht: hl.Table,
+) -> None:
+    mt = mt.select_cols(
+        discrepant=(
+            (
+                # All calls are diploid but the sex is Male
+                hl.agg.all(mt.GT.is_diploid())
+                & (sex_check_ht[mt.s].predicted_sex == Sex.MALE.value)
+            )
+            | (
+                # At least one call is haploid but the sex is Female
+                hl.agg.any(~mt.GT.is_diploid())
+                & (sex_check_ht[mt.s].predicted_sex == Sex.FEMALE.value)
+            )
+        ),
+    )
+    discrepant_rate = mt.aggregate_cols(hl.agg.fraction(mt.discrepant))
+    if discrepant_rate:
+        msg = f'{discrepant_rate:.2%} of samples have misaligned ploidy with their provided imputed sex.'
+        raise SeqrValidationError(msg)
+
+
 def validate_sample_type(
     mt: hl.MatrixTable,
     coding_and_noncoding_variants_ht: hl.Table,

v03_pipeline/lib/misc/validation_test.py

@@ -5,11 +5,14 @@
 from v03_pipeline.lib.misc.validation import (
     SeqrValidationError,
     validate_expected_contig_frequency,
+    validate_imputed_sex_ploidy,
     validate_no_duplicate_variants,
     validate_sample_type,
 )
 from v03_pipeline.lib.model import ReferenceGenome, SampleType
 
+TEST_SEX_CHECK_1 = 'v03_pipeline/var/test/sex_check/test_sex_check_1.ht'
+
 
 def _mt_from_contigs(contigs):
     return hl.MatrixTable.from_parts(
@@ -29,6 +32,57 @@ def _mt_from_contigs(contigs):
 
 
 class ValidationTest(unittest.TestCase):
+    def test_validate_imputed_sex_ploidy(self) -> None:
+        sex_check_ht = hl.read_table(TEST_SEX_CHECK_1)
+        mt = hl.MatrixTable.from_parts(
+            rows={
+                'locus': [
+                    hl.Locus(
+                        contig='chrX',
+                        position=1,
+                        reference_genome='GRCh38',
+                    ),
+                ],
+            },
+            cols={'s': ['HG00731_1', 'HG00732_1']},
+            entries={
+                'GT': [
+                    [
+                        hl.Call(alleles=[0, 0], phased=False),
+                        hl.Call(alleles=[0], phased=False),
+                    ],
+                ],
+            },
+        ).key_rows_by('locus')
+        validate_imputed_sex_ploidy(mt, sex_check_ht)
+        mt = hl.MatrixTable.from_parts(
+            rows={
+                'locus': [
+                    hl.Locus(
+                        contig='chrX',
+                        position=1,
+                        reference_genome='GRCh38',
+                    ),
+                ],
+            },
+            cols={'s': ['HG00731_1', 'HG00732_1']},
+            entries={
+                'GT': [
+                    [
+                        hl.Call(alleles=[0], phased=False),
+                        hl.Call(alleles=[0], phased=False),
+                    ],
+                ],
+            },
+        ).key_rows_by('locus')
+        self.assertRaisesRegex(
+            SeqrValidationError,
+            '50.00% of samples have misaligned ploidy',
+            validate_imputed_sex_ploidy,
+            mt,
+            sex_check_ht,
+        )
+
     def test_validate_no_duplicate_variants(self) -> None:
         mt = hl.MatrixTable.from_parts(
             rows={

v03_pipeline/lib/model/definitions.py

@@ -11,7 +11,6 @@ class AccessControl(Enum):
 class Sex(Enum):
     FEMALE = 'F'
     MALE = 'M'
-    UNKNOWN = 'U'
 
 
 class PipelineVersion(Enum):

v03_pipeline/lib/tasks/write_imported_callset.py

@@ -8,21 +8,25 @@
 )
 from v03_pipeline.lib.misc.validation import (
     validate_expected_contig_frequency,
+    validate_imputed_sex_ploidy,
     validate_no_duplicate_variants,
     validate_sample_type,
 )
 from v03_pipeline.lib.misc.vets import annotate_vets
 from v03_pipeline.lib.model import CachedReferenceDatasetQuery
 from v03_pipeline.lib.paths import (
     imported_callset_path,
+    sex_check_table_path,
     valid_cached_reference_dataset_query_path,
 )
 from v03_pipeline.lib.tasks.base.base_write import BaseWriteTask
 from v03_pipeline.lib.tasks.files import CallsetTask, GCSorLocalTarget, HailTableTask
+from v03_pipeline.lib.tasks.write_sex_check_table import WriteSexCheckTableTask
 
 
 class WriteImportedCallsetTask(BaseWriteTask):
     callset_path = luigi.Parameter()
+    imputed_sex_path = luigi.Parameter(default=None)
     filters_path = luigi.OptionalParameter(
         default=None,
         description='Optional path to part two outputs from callset (VCF shards containing filter information)',
@@ -35,6 +39,10 @@ class WriteImportedCallsetTask(BaseWriteTask):
         default=False,
         parsing=luigi.BoolParameter.EXPLICIT_PARSING,
     )
+    check_sex_and_relatedness = luigi.BoolParameter(
+        default=False,
+        parsing=luigi.BoolParameter.EXPLICIT_PARSING,
+    )
 
     def complete(self) -> luigi.Target:
         if not self.force and super().complete():
@@ -71,6 +79,20 @@ def requires(self) -> list[luigi.Task]:
                     ),
                 ),
             ]
+        if (
+            self.check_sex_and_relatedness
+            and self.dataset_type.check_sex_and_relatedness
+        ):
+            requirements = [
+                *requirements,
+                WriteSexCheckTableTask(
+                    self.reference_genome,
+                    self.dataset_type,
+                    self.sample_type,
+                    self.callset_path,
+                    self.imputed_sex_path,
+                ),
+            ]
         return [
             *requirements,
             CallsetTask(self.callset_path),
@@ -114,6 +136,21 @@ def create_table(self) -> hl.MatrixTable:
                 self.reference_genome,
                 self.sample_type,
             )
+        if (
+            self.check_sex_and_relatedness
+            and self.dataset_type.check_sex_and_relatedness
+        ):
+            sex_check_ht = hl.read_table(
+                sex_check_table_path(
+                    self.reference_genome,
+                    self.dataset_type,
+                    self.callset_path,
+                ),
+            )
+            validate_imputed_sex_ploidy(
+                mt,
+                sex_check_ht,
+            )
         return mt.annotate_globals(
             callset_path=self.callset_path,
             filters_path=self.filters_path or hl.missing(hl.tstr),

v03_pipeline/lib/tasks/write_remapped_and_subsetted_callset.py

@@ -65,16 +65,18 @@ def output(self) -> luigi.Target:
     def requires(self) -> list[luigi.Task]:
         requirements = [
             WriteImportedCallsetTask(
-                self.reference_genome,
-                self.dataset_type,
-                self.sample_type,
-                self.callset_path,
+                reference_genome=self.reference_genome,
+                dataset_type=self.dataset_type,
+                sample_type=self.sample_type,
+                callset_path=self.callset_path,
+                imputed_sex_path=self.imputed_sex_path,
                 # NB: filters_path is explicitly passed as None here
                 # to avoid carrying it throughout the rest of the pipeline.
                 # Only the primary import task itself should be aware of it.
-                None,
-                self.validate,
-                False,
+                filters_path=None,
+                validate=self.validate,
+                force=False,
+                check_sex_and_relatedness=self.check_sex_and_relatedness,
             ),
             RawFileTask(self.project_pedigree_path),
         ]

v03_pipeline/var/test/sex_check/test_imputed_sex_unexpected_value.tsv

@@ -0,0 +1,2 @@
+entity:sample_id	collaborator_participant_id	collaborator_sample_id	contamination_rate	coverage_region_1_metrics_file	crai_path	cram_md5_path	cram_path	datarepo_row_id	dragen_version	import:snapshot_id	import:timestamp	mapped_percentage	mapping_metrics_file	material_type	mean_coverage	original_material_type	participant_id	pass_fail_value	pdo	percent_bases_at_20x	percent_callability	predicted_sex	product	receipt_date	reported_sex	research_project	single_sample_vcf_index_path	single_sample_vcf_md5_path	single_sample_vcf_path	total_bases	variant_calling_metrics_file
+SM-DM69X	abc_2	abc_2	0E+00	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/c4c07edf-7735-4aa7-9283-7cb2607b60a2/GLE-5774-3-3.qc-coverage-region-1_coverage_metrics.csv	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/dcd4c271-0249-47f1-8e91-81f74735c5a1/GLE-5774-3-3.cram.crai	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/ec41ec06-673f-4fe2-a063-23dc5fe1dcce/GLE-5774-3-3.cram.md5sum	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/aad0e270-2ad5-4f39-b968-9b4beafeb5cc/GLE-5774-3-3.cram	a4b04a39-9234-4028-a155-442c4acf12a0	07.021.604.3.7.8	ce74d94c-c33d-49d7-85c9-5f3cbd08aff7	2024-04-17T15:02:46	99.800000000	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/c3a9e6f2-4c68-410b-823d-46ca406e5061/GLE-5774-3-3.mapping_metrics.csv	DNA:DNA Genomic	35.300000000	Whole Blood:Whole Blood	PT-24OHM	Pass	PDO-32755	96.320000000	97.340000000	Unknown	P-WG-0139	2017-04-12 04:00:00	Unknown	RP-3061	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/c71cd2a1-c789-4715-9ebc-dbfc40d9f2e2/GLE-5774-3-3.vcf.gz.tbi	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/957a99cb-c9a9-4fc5-a0ec-53f9e461469e/GLE-5774-3-3.vcf.gz.md5sum	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/df520949-5f2b-4976-9d46-80d1cc299813/GLE-5774-3-3.vcf.gz	133253714921.000000000	gs://datarepo-556a9c15-bucket/2a4202b0-93f5-4ebe-8d2b-fd4cfb2b881d/2e98e51b-9394-4e64-977f-e9010a4e16dc/GLE-5774-3-3.vc_metrics.csv