merge

Author: bpblanken

v03_pipeline/lib/misc/terra_data_repository.py

@@ -31,7 +31,7 @@ def _tdr_request(resource: str) -> dict:
 
 
 def _get_dataset_ids() -> list[str]:
-    res_body = _tdr_request('datasets')
+    res_body = _tdr_request('datasets?limit=50000')  # Arbitrary large number :/
     items = res_body['items']
     for item in items:
         if not any(x['cloudResource'] == BIGQUERY_RESOURCE for x in item['storage']):

v03_pipeline/lib/misc/terra_data_repository_test.py

@@ -137,7 +137,7 @@ class TerraDataRepositoryTest(unittest.TestCase):
     @responses.activate
     def test_get_dataset_ids(self, _: Mock) -> None:
         responses.get(
-            os.path.join(TDR_ROOT_URL, 'datasets'),
+            os.path.join(TDR_ROOT_URL, 'datasets?limit=50000'),
             body=json.dumps(
                 {
                     'total': 3,

v03_pipeline/lib/reference_datasets/exac_test.py

@@ -52,3 +52,9 @@ def test_exac(self):
                     ),
                 ],
             )
+            self.assertEqual(
+                ht.globals.collect(),
+                [
+                    hl.Struct(version='1.1', enums=hl.Struct()),
+                ],
+            )

v03_pipeline/lib/reference_datasets/gnomad_mito.py

@@ -9,6 +9,6 @@ def get_ht(path: str, *_) -> hl.Table:
         AC_het=hl.int32(ht.AC_het),
         AF_hom=ht.AF_hom,
         AF_het=ht.AF_het,
-        max_hl=ht.max_hl,
+        max_hl=hl.float32(ht.max_hl),
     )
     return ht.select_globals()

v03_pipeline/lib/reference_datasets/reference_dataset.py

@@ -90,7 +90,7 @@ def enums(self) -> dict | None:
     def enum_globals(self) -> hl.Struct:
         if self.enums:
             return hl.Struct(**self.enums)
-        return hl.missing(hl.tstruct(hl.tstr, hl.tarray(hl.tstr)))
+        return hl.Struct()
 
     @property
     def filter(  # noqa: A003
@@ -263,15 +263,15 @@ def get_ht(
         },
         ReferenceGenome.GRCh37: {
             DATASET_TYPES: frozenset([DatasetType.SNV_INDEL]),
-            VERSION: '1.0',
+            VERSION: '1.1',
             PATH: [
                 'gs://seqr-reference-data/GRCh37/spliceai/spliceai_scores.masked.snv.hg19.vcf.gz',
                 'gs://seqr-reference-data/GRCh37/spliceai/spliceai_scores.masked.indel.hg19.vcf.gz',
             ],
         },
         ReferenceGenome.GRCh38: {
             DATASET_TYPES: frozenset([DatasetType.SNV_INDEL]),
-            VERSION: '1.0',
+            VERSION: '1.1',
             # NB: SpliceAI data is only available to download for authenticated Illumina users, so we will host the data
             PATH: [
                 'gs://seqr-reference-data/GRCh38/spliceai/spliceai_scores.masked.snv.hg38.vcf.gz',
@@ -369,7 +369,7 @@ def get_ht(
     ReferenceDataset.gnomad_mito: {
         ReferenceGenome.GRCh38: {
             DATASET_TYPES: frozenset([DatasetType.MITO]),
-            VERSION: '1.0',
+            VERSION: '1.1',
             PATH: 'gs://gcp-public-data--gnomad/release/3.1/ht/genomes/gnomad.genomes.v3.1.sites.chrM.ht',
         },
     },

v03_pipeline/lib/reference_datasets/splice_ai.py

@@ -9,6 +9,26 @@
 from v03_pipeline.lib.reference_datasets.misc import vcf_to_ht
 
 
+def remove_duplicate_scores(ht: hl.Table):
+    #
+    # SpliceAI has many duplicate rows of the ilk:
+    #
+    #  1:861264      | ["C","A"]  | NA   | -1.00e+01 | NA       | ["A|AL645608.1|0.00|0.00|0.00|0.00|2|27|12|1"] |
+    #  1:861264      | ["C","A"]  | NA   | -1.00e+01 | NA       | ["A|SAMD11|0.02|0.01|0.00|0.00|14|38|14|38"]
+    #
+    count_ht = ht.group_by(*ht.key).aggregate(n=hl.agg.count())
+    duplicate_variants_ht = count_ht.filter(count_ht.n > 1)
+    duplicates_ht = ht.semi_join(duplicate_variants_ht)
+    non_duplicates_ht = ht.anti_join(duplicates_ht)
+    return non_duplicates_ht.union(
+        # Remove rows that 1) are part of duplicate variant groupings
+        # and 2) contain dots.  Then, remove arbitrarily with .distinct()
+        duplicates_ht.filter(
+            ~duplicates_ht.info.SpliceAI[0].split(delim='\\|')[1].contains('.'),
+        ).distinct(),
+    )
+
+
 def get_ht(
     paths: list[str],
     reference_genome: ReferenceGenome,
@@ -26,6 +46,7 @@ def get_ht(
         # of partititons.
     )
     ht, _ = checkpoint(ht)
+    ht = remove_duplicate_scores(ht)
 
     # SpliceAI INFO field description from the VCF header: SpliceAIv1.3 variant annotation. These include
     # delta scores (DS) and delta positions (DP) for acceptor gain (AG), acceptor loss (AL), donor gain (DG), and
@@ -39,7 +60,7 @@ def get_ht(
         .map(hl.float32),
     )
     ht = ht.annotate(delta_score=hl.max(ht.delta_scores))
-    ht = ht.annotate(
+    return ht.annotate(
         splice_consequence_id=hl.if_else(
             ht.delta_score > 0,
             # Splice Consequence enum ID is the index of the max score
@@ -48,6 +69,3 @@ def get_ht(
             num_delta_scores,
         ),
     ).drop('delta_scores')
-    return ht.group_by(*ht.key).aggregate(
-        splice_consequence_id=hl.agg.min(ht.splice_consequence_id),
-    )

v03_pipeline/lib/reference_datasets/topmed.py

@@ -15,7 +15,8 @@
 
 def get_ht(path: str, reference_genome: ReferenceGenome) -> hl.Table:
     ht = vcf_to_ht(path, reference_genome)
-    ht = ht.filter(ht.locus.position == hl.int(ht.info.OriginalStart))
+    if reference_genome == ReferenceGenome.GRCh37:
+        ht = ht.filter(ht.locus.position == hl.int(ht.info.OriginalStart))
     return ht.select(
         **{k: parse_nested_field(ht, v) for k, v in SELECT.items()},
     )

v03_pipeline/lib/tasks/reference_data/update_variant_annotations_table_with_updated_reference_dataset_test.py

@@ -146,7 +146,7 @@ def test_update_vat_snv_indel_38(
                             eigen='1.1',
                             clinvar='2024-11-11',
                             exac='1.1',
-                            splice_ai='1.0',
+                            splice_ai='1.1',
                             topmed='1.1',
                             hgmd='1.0',
                             gnomad_exomes='1.0',
@@ -297,7 +297,7 @@ def test_update_vat_mito_38(
                             hmtvar='1.1',
                             mitimpact='1.0',
                             mitomap='1.0',
-                            gnomad_mito='1.0',
+                            gnomad_mito='1.1',
                             local_constraint_mito='1.0',
                             clinvar='2024-11-11',
                             dbnsfp='1.0',
@@ -425,7 +425,7 @@ def test_update_vat_snv_indel_37(
                             eigen='1.1',
                             clinvar='2024-11-11',
                             exac='1.1',
-                            splice_ai='1.0',
+                            splice_ai='1.1',
                             topmed='1.1',
                             hgmd='1.0',
                             gnomad_exomes='1.0',

v03_pipeline/lib/tasks/update_variant_annotations_table_with_new_samples_test.py

@@ -507,7 +507,7 @@ def test_multiple_update_vat(
                             exac='1.1',
                             gnomad_exomes='1.0',
                             gnomad_genomes='1.0',
-                            splice_ai='1.0',
+                            splice_ai='1.1',
                             topmed='1.1',
                             gnomad_non_coding_constraint='1.0',
                             screen='1.0',
@@ -772,7 +772,7 @@ def test_update_vat_without_accessing_private_datasets(
                         exac='1.1',
                         gnomad_exomes='1.0',
                         gnomad_genomes='1.0',
-                        splice_ai='1.0',
+                        splice_ai='1.1',
                         topmed='1.1',
                         gnomad_non_coding_constraint='1.0',
                         screen='1.0',
@@ -827,7 +827,7 @@ def test_mito_update_vat(
                         versions=hl.Struct(
                             clinvar='2024-11-11',
                             dbnsfp='1.0',
-                            gnomad_mito='1.0',
+                            gnomad_mito='1.1',
                             helix_mito='1.0',
                             hmtvar='1.1',
                             mitomap='1.0',

v03_pipeline/lib/tasks/write_tdr_metrics_file.py

@@ -1,8 +1,10 @@
 import csv
 
+import google.api_core.exceptions
 import luigi
 import luigi.util
 
+from v03_pipeline.lib.logger import get_logger
 from v03_pipeline.lib.misc.terra_data_repository import (
     BIGQUERY_METRICS,
     bq_metrics_query,
@@ -13,6 +15,8 @@
 )
 from v03_pipeline.lib.tasks.files import GCSorLocalTarget
 
+logger = get_logger(__name__)
+
 
 @luigi.util.inherits(BaseLoadingPipelineParams)
 class WriteTDRMetricsFileTask(luigi.Task):
@@ -31,5 +35,8 @@ def run(self):
         with self.output().open('w') as f:
             writer = csv.DictWriter(f, fieldnames=BIGQUERY_METRICS, delimiter='\t')
             writer.writeheader()
-            for row in bq_metrics_query(self.bq_table_name):
-                writer.writerow(row)
+            try:
+                for row in bq_metrics_query(self.bq_table_name):
+                    writer.writerow(row)
+            except google.api_core.exceptions.BadRequest:
+                logger.exception('Query Failed')