Merge pull request #779 from broadinstitute/dev

main <- dev

Author: bpblanken

v03_pipeline/bin/vep-110-GRCh38.sh

@@ -12,7 +12,7 @@ export PROJECT="$(gcloud config get-value project)"
 export VEP_CONFIG_PATH="$(/usr/share/google/get_metadata_value attributes/VEP_CONFIG_PATH)"
 export VEP_REPLICATE="$(/usr/share/google/get_metadata_value attributes/VEP_REPLICATE)"
 export ASSEMBLY=GRCh38
-export VEP_DOCKER_IMAGE=gcr.io/seqr-project/vep-docker-image
+export VEP_DOCKER_IMAGE=gcr.io/seqr-project/vep-docker-image:110
 
 mkdir -p /vep_data
 
@@ -61,7 +61,7 @@ gcloud storage cat --billing-project $PROJECT gs://seqr-reference-data/vep/110/h
 # gzip -d Homo_sapiens.${ASSEMBLY}.dna.primary_assembly.fa.gz
 # bgzip Homo_sapiens.${ASSEMBLY}.dna.primary_assembly.fa
 # samtools faidx Homo_sapiens.${ASSEMBLY}.dna.primary_assembly.fa.gz
-gcloud storage cp --billing-project $PROJECT 'gs://seqr-reference-data/vep/110/Homo_sapiens.${ASSEMBLY}.dna.primary_assembly.fa.*' /vep_data/ &
+gcloud storage cp --billing-project $PROJECT "gs://seqr-reference-data/vep/110/Homo_sapiens.${ASSEMBLY}.dna.primary_assembly.fa.*" /vep_data/ &
 docker pull ${VEP_DOCKER_IMAGE} &
 wait
 
@@ -88,3 +88,4 @@ docker run -i -v /vep_data/:/opt/vep/.vep/:ro ${VEP_DOCKER_IMAGE} \
   /opt/vep/src/ensembl-vep/vep "\$@"
 EOF
 chmod +x /vep.sh
+

v03_pipeline/lib/methods/sex_check.py

@@ -6,13 +6,14 @@
 from v03_pipeline.lib.paths import remapped_and_subsetted_callset_path
 
 
-def get_callset_ht(
+def get_callset_ht(  # noqa: PLR0913
     reference_genome: ReferenceGenome,
     dataset_type: DatasetType,
     callset_paths: list[str],
     project_guids: list[str],
     project_remap_paths: list[str],
     project_pedigree_paths: list[str],
+    imputed_sex_paths: list[str] | None,
 ):
     callset_hts = [
         hl.read_matrix_table(
@@ -23,11 +24,12 @@ def get_callset_ht(
                 project_guid,
             ),
         ).rows()
-        for (callset_path, project_guid, _, _) in callset_project_pairs(
+        for (callset_path, project_guid, _, _, _) in callset_project_pairs(
             callset_paths,
             project_guids,
             project_remap_paths,
             project_pedigree_paths,
+            imputed_sex_paths,
         )
     ]
 
@@ -49,18 +51,34 @@ def callset_project_pairs(
     project_guids: list[str],
     project_remap_paths: list[str],
     project_pedigree_paths: list[str],
+    imputed_sex_paths: list[str] | None,
 ):
     if len(callset_paths) == len(project_guids):
         return zip(
             callset_paths,
             project_guids,
             project_remap_paths,
             project_pedigree_paths,
+            imputed_sex_paths
+            if imputed_sex_paths is not None
+            else [None] * len(callset_paths),
             strict=True,
         )
     return (
-        (callset_path, project_guid, project_remap_path, project_pedigree_path)
-        for callset_path in callset_paths
+        (
+            callset_path,
+            project_guid,
+            project_remap_path,
+            project_pedigree_path,
+            imputed_sex_path,
+        )
+        for callset_path, imputed_sex_path in zip(
+            callset_paths,
+            imputed_sex_paths
+            if imputed_sex_paths is not None
+            else [None] * len(callset_paths),
+            strict=False,
+        )
         for (project_guid, project_remap_path, project_pedigree_path) in zip(
             project_guids,
             project_remap_paths,

v03_pipeline/lib/methods/sex_check_test.py

@@ -107,7 +107,6 @@ def build_relatedness_check_lookup(
     relatedness_check_ht: hl.Table,
     remap_lookup: hl.dict,
 ) -> dict[tuple[str, str], list]:
-    # Build relatedness check lookup
     relatedness_check_ht = relatedness_check_ht.key_by(
         i=remap_lookup.get(relatedness_check_ht.i, relatedness_check_ht.i),
         j=remap_lookup.get(relatedness_check_ht.j, relatedness_check_ht.j),
@@ -124,12 +123,11 @@ def build_sex_check_lookup(
     sex_check_ht: hl.Table,
     remap_lookup: hl.dict,
 ) -> dict[str, Sex]:
-    # Build sex check lookup
     sex_check_ht = sex_check_ht.key_by(
         s=remap_lookup.get(sex_check_ht.s, sex_check_ht.s),
     )
-    sex_check_ht = sex_check_ht.select('sex')
-    return {r.s: Sex(r.sex) for r in sex_check_ht.collect()}
+    sex_check_ht = sex_check_ht.select('predicted_sex')
+    return {r.s: Sex(r.predicted_sex) for r in sex_check_ht.collect()}
 
 
 def get_families_failed_missing_samples(

v03_pipeline/lib/misc/callsets.py

@@ -56,16 +56,16 @@ def test_build_relatedness_check_lookup(self):
     def test_build_sex_check_lookup(self):
         ht = hl.Table.parallelize(
             [
-                {'s': 'remapped_id', 'sex': 'M'},
-                {'s': 'ROS_006_18Y03227_D1', 'sex': 'M'},
-                {'s': 'ROS_006_18Y03228_D1', 'sex': 'M'},
-                {'s': 'ROS_007_19Y05919_D1', 'sex': 'M'},
-                {'s': 'ROS_007_19Y05939_D1', 'sex': 'F'},
-                {'s': 'ROS_007_19Y05987_D1', 'sex': 'M'},
+                {'s': 'remapped_id', 'predicted_sex': 'M'},
+                {'s': 'ROS_006_18Y03227_D1', 'predicted_sex': 'M'},
+                {'s': 'ROS_006_18Y03228_D1', 'predicted_sex': 'M'},
+                {'s': 'ROS_007_19Y05919_D1', 'predicted_sex': 'M'},
+                {'s': 'ROS_007_19Y05939_D1', 'predicted_sex': 'F'},
+                {'s': 'ROS_007_19Y05987_D1', 'predicted_sex': 'M'},
             ],
             hl.tstruct(
                 s=hl.tstr,
-                sex=hl.tstr,
+                predicted_sex=hl.tstr,
             ),
             key='s',
         )
@@ -178,16 +178,16 @@ def test_all_relatedness_checks(self):
     def test_get_families_failed_sex_check(self):
         sex_check_ht = hl.Table.parallelize(
             [
-                {'s': 'ROS_006_18Y03226_D1', 'sex': 'M'},
-                {'s': 'ROS_006_18Y03227_D1', 'sex': 'F'},
-                {'s': 'ROS_006_18Y03228_D1', 'sex': 'F'},
-                {'s': 'ROS_007_19Y05919_D1', 'sex': 'F'},
-                {'s': 'ROS_007_19Y05939_D1', 'sex': 'F'},
-                {'s': 'ROS_007_19Y05987_D1', 'sex': 'F'},
+                {'s': 'ROS_006_18Y03226_D1', 'predicted_sex': 'M'},
+                {'s': 'ROS_006_18Y03227_D1', 'predicted_sex': 'F'},
+                {'s': 'ROS_006_18Y03228_D1', 'predicted_sex': 'F'},
+                {'s': 'ROS_007_19Y05919_D1', 'predicted_sex': 'F'},
+                {'s': 'ROS_007_19Y05939_D1', 'predicted_sex': 'F'},
+                {'s': 'ROS_007_19Y05987_D1', 'predicted_sex': 'F'},
             ],
             hl.tstruct(
                 s=hl.tstr,
-                sex=hl.tstr,
+                predicted_sex=hl.tstr,
             ),
             key='s',
         )

v03_pipeline/lib/misc/family_loading_failures.py

@@ -5,12 +5,15 @@
 import hail as hl
 
 from v03_pipeline.lib.misc.gcnv import parse_gcnv_genes
-from v03_pipeline.lib.model import DatasetType, Env, ReferenceGenome
+from v03_pipeline.lib.model import DatasetType, Env, ReferenceGenome, Sex
 
 BIALLELIC = 2
 B_PER_MB = 1 << 20  # 1024 * 1024
 MB_PER_PARTITION = 128
 
+MALE = 'Male'
+FEMALE = 'Female'
+
 
 def does_file_exist(path: str) -> bool:
     if path.startswith('gs://'):
@@ -156,6 +159,25 @@ def select_relevant_fields(
     return mt.select_entries(*dataset_type.entries_fields)
 
 
+def import_imputed_sex(imputed_sex_path: str) -> hl.Table:
+    ht = hl.import_table(imputed_sex_path)
+    ht = ht.select(
+        s=ht.collaborator_sample_id,
+        predicted_sex=(
+            hl.case()
+            .when(ht.predicted_sex == FEMALE, Sex.FEMALE.value)
+            .when(ht.predicted_sex == MALE, Sex.MALE.value)
+            .or_error(
+                hl.format(
+                    'Found unexpected value %s in imputed sex file',
+                    ht.predicted_sex,
+                ),
+            )
+        ),
+    )
+    return ht.key_by(ht.s)
+
+
 def import_remap(remap_path: str) -> hl.Table:
     ht = hl.import_table(remap_path)
     ht = ht.select(

v03_pipeline/lib/misc/family_loading_failures_test.py

@@ -1,7 +1,17 @@
 import unittest
 
-from v03_pipeline.lib.misc.io import compute_hail_n_partitions, file_size_bytes
+import hail as hl
 
+from v03_pipeline.lib.misc.io import (
+    compute_hail_n_partitions,
+    file_size_bytes,
+    import_imputed_sex,
+)
+
+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'
 
@@ -17,3 +27,22 @@ def test_compute_hail_n_partitions(self) -> None:
         self.assertEqual(compute_hail_n_partitions(23), 1)
         self.assertEqual(compute_hail_n_partitions(191310), 1)
         self.assertEqual(compute_hail_n_partitions(1913100000), 15)
+
+    def test_import_imputed_sex(self) -> None:
+        ht = import_imputed_sex(TEST_IMPUTED_SEX)
+        self.assertListEqual(
+            ht.collect(),
+            [
+                hl.Struct(s='abc_1', predicted_sex='M'),
+                hl.Struct(s='abc_2', predicted_sex='F'),
+                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,
+        )