feat: add gnomad_svs as reference dataset. (#1041)

* bump hail version (#1027)

* First pass

* Finish test

* do it correctly :/

* ruff

* correct feature flag

* bump

* fix test

* fix svs

* ruff

* Benb/bugfix allow deleting project to succeed if no project (#1038)

* First pass

* add comment

* ruff

* gnomad svs reference dataset

* tweaks

* tdr updates for sample qc

* dead code (#1045)

* add global to table

* ruff

* use the function

* add missing enums struct

* sprawl :/

* test finished!

* ruff

* change field name

* delete

* fixes

* need loadable families

* ruff

---------

Co-authored-by: Julia Klugherz <juliaklugherz@gmail.com>

Author: bpblanken

v03_pipeline/lib/annotations/fields_test.py

@@ -131,7 +131,7 @@ def test_get_formatting_fields(self, mock_vep: Mock) -> None:
                                 reference_genome,
                                 DatasetType.SNV_INDEL,
                             )
-                            if reference_dataset.is_keyed_by_interval
+                            if reference_dataset.formatting_annotation
                         },
                         **(
                             {

v03_pipeline/lib/annotations/sv.py

@@ -175,11 +175,14 @@ def end_locus(ht: hl.Table, **_: Any) -> hl.StructExpression:
     )
 
 
-def gnomad_svs(ht: hl.Table, **_: Any) -> hl.Expression:
-    return hl.or_missing(
-        hl.is_defined(ht['info.gnomAD_V2_AF']),
-        hl.struct(AF=hl.float32(ht['info.gnomAD_V2_AF']), ID=ht['info.gnomAD_V2_SVID']),
-    )
+def gnomad_svs(
+    ht: hl.Table,
+    gnomad_svs_ht: hl.Table,
+    **_: Any,
+) -> hl.Expression:
+    return gnomad_svs_ht.annotate(
+        ID=gnomad_svs_ht.KEY,
+    )[ht['info.GNOMAD_V4.1_TRUTH_VID']]
 
 
 def gt_stats(ht: hl.Table, **_: Any) -> hl.Expression:

v03_pipeline/lib/misc/io_test.py

@@ -23,7 +23,6 @@
 TEST_PEDIGREE_3 = 'v03_pipeline/var/test/pedigrees/test_pedigree_3.tsv'
 TEST_MITO_MT = 'v03_pipeline/var/test/callsets/mito_1.mt'
 TEST_REMAP = 'v03_pipeline/var/test/remaps/test_remap_1.tsv'
-TEST_SV_VCF = 'v03_pipeline/var/test/callsets/sv_1.vcf'
 
 
 class IOTest(unittest.TestCase):

v03_pipeline/lib/model/dataset_type.py

@@ -121,10 +121,9 @@ def row_fields(
                 'info.CPX_TYPE': hl.tstr,
                 'info.END': hl.tint32,
                 'info.END2': hl.tint32,
-                'info.gnomAD_V2_AF': hl.tfloat64,
-                'info.gnomAD_V2_SVID': hl.tstr,
                 'info.N_HET': hl.tint32,
                 'info.N_HOMALT': hl.tint32,
+                'info.GNOMAD_V4.1_TRUTH_VID': hl.tstr,
                 'info.StrVCTVRE': hl.tstr,
                 'info.SVLEN': hl.tint32,
                 **sv.CONSEQ_PREDICTED_GENE_COLS,

v03_pipeline/lib/reference_datasets/gnomad_svs.py

@@ -0,0 +1,18 @@
+import hail as hl
+
+from v03_pipeline.lib.model import ReferenceGenome
+from v03_pipeline.lib.reference_datasets.misc import vcf_to_ht
+
+
+def get_ht(path: str, reference_genome: ReferenceGenome) -> hl.Table:
+    ht = vcf_to_ht(path, reference_genome)
+    ht = ht.select(
+        KEY=ht.rsid,
+        AF=ht.info.AF[0],
+        AC=ht.info.AC[0],
+        AN=ht.info.AN,
+        N_HET=ht.info.N_HET,
+        N_HOMREF=ht.info.N_HOMREF,
+    )
+    ht = ht.key_by('KEY')
+    return ht.drop('locus', 'alleles')

v03_pipeline/lib/reference_datasets/misc.py

@@ -83,9 +83,12 @@ def filter_contigs(ht, reference_genome: ReferenceGenome):
                 ht.interval.start.contig,
             ),
         )
-    return ht.filter(
-        hl.set(reference_genome.standard_contigs).contains(ht.locus.contig),
-    )
+    # SV reference datasets are not keyed by locus.
+    if hasattr(ht, 'locus'):
+        return ht.filter(
+            hl.set(reference_genome.standard_contigs).contains(ht.locus.contig),
+        )
+    return ht
 
 
 def vcf_to_ht(

v03_pipeline/lib/reference_datasets/reference_dataset.py

@@ -6,6 +6,7 @@
 
 import hail as hl
 
+from v03_pipeline.lib.annotations import snv_indel, sv
 from v03_pipeline.lib.misc.validation import (
     validate_allele_type,
     validate_no_duplicate_variants,
@@ -27,8 +28,9 @@
 DATASET_TYPES = 'dataset_types'
 ENUMS = 'enums'
 EXCLUDE_FROM_ANNOTATIONS = 'exclude_from_annotations'
+EXCLUDE_FROM_ANNOTATIONS_UPDATES = 'exclude_from_annotations_updates'
+FORMATTING_ANNOTATION = 'formatting_annotation'
 FILTER = 'filter'
-IS_INTERVAL = 'is_interval'
 SELECT = 'select'
 VERSION = 'version'
 PATH = 'path'
@@ -69,9 +71,24 @@ def for_reference_genome_dataset_type_annotations(
             if not CONFIG[dataset].get(EXCLUDE_FROM_ANNOTATIONS, False)
         }
 
+    @classmethod
+    def for_reference_genome_dataset_type_annotations_updates(
+        cls,
+        reference_genome: ReferenceGenome,
+        dataset_type: DatasetType,
+    ) -> set['ReferenceDataset']:
+        return {
+            dataset
+            for dataset in cls.for_reference_genome_dataset_type_annotations(
+                reference_genome,
+                dataset_type,
+            )
+            if not CONFIG[dataset].get(EXCLUDE_FROM_ANNOTATIONS_UPDATES, False)
+        }
+
     @property
-    def is_keyed_by_interval(self) -> bool:
-        return CONFIG[self].get(IS_INTERVAL, False)
+    def formatting_annotation(self) -> Callable | None:
+        return CONFIG[self].get(FORMATTING_ANNOTATION)
 
     @property
     def access_control(self) -> AccessControl:
@@ -154,9 +171,10 @@ class ReferenceDataset(BaseReferenceDataset, StrEnum):
     gnomad_coding_and_noncoding = 'gnomad_coding_and_noncoding'
     gnomad_exomes = 'gnomad_exomes'
     gnomad_genomes = 'gnomad_genomes'
-    gnomad_qc = 'gnomad_qc'
     gnomad_mito = 'gnomad_mito'
     gnomad_non_coding_constraint = 'gnomad_non_coding_constraint'
+    gnomad_qc = 'gnomad_qc'
+    gnomad_svs = 'gnomad_svs'
     screen = 'screen'
     local_constraint_mito = 'local_constraint_mito'
     mitomap = 'mitomap'
@@ -394,7 +412,7 @@ def get_ht(
         },
     },
     ReferenceDataset.gnomad_non_coding_constraint: {
-        IS_INTERVAL: True,
+        FORMATTING_ANNOTATION: snv_indel.gnomad_non_coding_constraint,
         ReferenceGenome.GRCh38: {
             DATASET_TYPES: frozenset([DatasetType.SNV_INDEL]),
             VERSION: '1.0',
@@ -414,7 +432,7 @@ def get_ht(
                 'low-DNase',
             ],
         },
-        IS_INTERVAL: True,
+        FORMATTING_ANNOTATION: snv_indel.screen,
         ReferenceGenome.GRCh38: {
             DATASET_TYPES: frozenset([DatasetType.SNV_INDEL]),
             VERSION: '1.0',
@@ -428,6 +446,15 @@ def get_ht(
             PATH: 'https://www.biorxiv.org/content/biorxiv/early/2023/01/27/2022.12.16.520778/DC3/embed/media-3.zip',
         },
     },
+    ReferenceDataset.gnomad_svs: {
+        EXCLUDE_FROM_ANNOTATIONS_UPDATES: True,
+        FORMATTING_ANNOTATION: sv.gnomad_svs,
+        ReferenceGenome.GRCh38: {
+            DATASET_TYPES: frozenset([DatasetType.SV]),
+            VERSION: '1.0',
+            PATH: 'gs://gcp-public-data--gnomad/release/4.1/genome_sv/gnomad.v4.1.sv.sites.vcf.gz',
+        },
+    },
 }
 CONFIG[ReferenceDatasetQuery.clinvar_path_variants] = {
     EXCLUDE_FROM_ANNOTATIONS: True,

v03_pipeline/lib/tasks/base/base_update_variant_annotations_table.py

@@ -8,6 +8,7 @@
 )
 from v03_pipeline.lib.reference_datasets.reference_dataset import (
     BaseReferenceDataset,
+    ReferenceDataset,
     ReferenceDatasetQuery,
 )
 from v03_pipeline.lib.tasks.base.base_update import BaseUpdateTask
@@ -77,17 +78,13 @@ def update_table(self, ht: hl.Table) -> hl.Table:
     def annotate_globals(
         self,
         ht: hl.Table,
+        reference_datasets: set[ReferenceDataset],
     ) -> hl.Table:
         ht = ht.annotate_globals(
             versions=hl.Struct(),
             enums=hl.Struct(),
         )
-        for (
-            reference_dataset
-        ) in BaseReferenceDataset.for_reference_genome_dataset_type_annotations(
-            self.reference_genome,
-            self.dataset_type,
-        ):
+        for reference_dataset in reference_datasets:
             rd_ht = hl.read_table(
                 valid_reference_dataset_path(self.reference_genome, reference_dataset),
             )

v03_pipeline/lib/tasks/reference_data/update_variant_annotations_table_with_updated_reference_dataset.py

@@ -29,7 +29,7 @@ def __init__(self, *args, **kwargs):
     def complete(self) -> bool:
         reference_dataset_names = {
             rd.name
-            for rd in BaseReferenceDataset.for_reference_genome_dataset_type_annotations(
+            for rd in BaseReferenceDataset.for_reference_genome_dataset_type_annotations_updates(
                 self.reference_genome,
                 self.dataset_type,
             )
@@ -65,18 +65,11 @@ def update_table(self, ht: hl.Table) -> hl.Table:
             reference_dataset_ht = hl.read_table(
                 valid_reference_dataset_path(self.reference_genome, reference_dataset),
             )
-            if reference_dataset.is_keyed_by_interval:
-                formatting_fn = next(
-                    x
-                    for x in self.dataset_type.formatting_annotation_fns(
-                        self.reference_genome,
-                    )
-                    if x.__name__ == reference_dataset.name
-                )
+            if reference_dataset.formatting_annotation:
                 ht = ht.annotate(
                     **get_fields(
                         ht,
-                        [formatting_fn],
+                        [reference_dataset.formatting_annotation],
                         **{f'{reference_dataset.name}_ht': reference_dataset_ht},
                         **self.param_kwargs,
                     ),
@@ -103,4 +96,10 @@ def update_table(self, ht: hl.Table) -> hl.Table:
                 )
                 ht = ht.join(reference_dataset_ht, 'left')
 
-        return self.annotate_globals(ht)
+        return self.annotate_globals(
+            ht,
+            BaseReferenceDataset.for_reference_genome_dataset_type_annotations_updates(
+                self.reference_genome,
+                self.dataset_type,
+            ),
+        )

v03_pipeline/lib/tasks/update_variant_annotations_table_with_new_samples.py

@@ -9,6 +9,9 @@
     lookup_table_path,
     new_variants_table_path,
 )
+from v03_pipeline.lib.reference_datasets.reference_dataset import (
+    BaseReferenceDataset,
+)
 from v03_pipeline.lib.tasks.base.base_loading_run_params import (
     BaseLoadingRunParams,
 )
@@ -104,7 +107,13 @@ def update_table(self, ht: hl.Table) -> hl.Table:
             ht = non_callset_variants_ht.union(callset_variants_ht, unify=True)
 
         # Fix up the globals and mark the table as updated with these callset/project pairs.
-        ht = self.annotate_globals(ht)
+        ht = self.annotate_globals(
+            ht,
+            BaseReferenceDataset.for_reference_genome_dataset_type_annotations(
+                self.reference_genome,
+                self.dataset_type,
+            ),
+        )
         return ht.annotate_globals(
             updates=ht.updates.union(
                 {

v03_pipeline/lib/tasks/update_variant_annotations_table_with_new_samples_test.py

@@ -935,8 +935,9 @@ def test_sv_update_vat(
             ht.globals.collect(),
             [
                 hl.Struct(
-                    versions=hl.Struct(),
+                    versions=hl.Struct(gnomad_svs='1.0'),
                     enums=hl.Struct(
+                        gnomad_svs=hl.Struct(),
                         sv_type=SV_TYPES,
                         sv_type_detail=SV_TYPE_DETAILS,
                         sorted_gene_consequences=hl.Struct(
@@ -1548,8 +1549,12 @@ def test_sv_update_vat(
                         Het=590,
                     ),
                     gnomad_svs=hl.Struct(
-                        AF=hl.eval(hl.float32(0.06896299868822098)),
-                        ID='gnomAD-SV_v2.1_INS_chr1_65',
+                        ID='gnomAD-SV_v3_BND_chr1_0001f5b7',
+                        AF=3.2e-05,
+                        AC=4,
+                        AN=126092,
+                        N_HET=4,
+                        N_HOMREF=63042,
                     ),
                     rg37_locus=hl.Locus(
                         contig=1,

v03_pipeline/lib/tasks/write_new_variants_table.py

@@ -168,7 +168,7 @@ def create_table(self) -> hl.Table:
             self.reference_genome,
             self.dataset_type,
         ):
-            if reference_dataset.is_keyed_by_interval:
+            if reference_dataset.formatting_annotation:
                 continue
             reference_dataset_ht = self.annotation_dependencies[
                 f'{reference_dataset.value}_ht'

v03_pipeline/lib/tasks/write_variant_annotations_vcf_test.py

@@ -10,8 +10,9 @@
 from v03_pipeline.lib.tasks.write_variant_annotations_vcf import (
     WriteVariantAnnotationsVCF,
 )
-from v03_pipeline.lib.test.mock_complete_task import MockCompleteTask
-from v03_pipeline.lib.test.mocked_dataroot_testcase import MockedDatarootTestCase
+from v03_pipeline.lib.test.mocked_reference_datasets_testcase import (
+    MockedReferenceDatasetsTestCase,
+)
 
 TEST_SV_VCF = 'v03_pipeline/var/test/callsets/sv_1.vcf'
 TEST_PEDIGREE_5 = 'v03_pipeline/var/test/pedigrees/test_pedigree_5.tsv'
@@ -36,20 +37,15 @@
 }
 
 
-class WriteVariantAnnotationsVCFTest(MockedDatarootTestCase):
+class WriteVariantAnnotationsVCFTest(MockedReferenceDatasetsTestCase):
     @patch(
         'v03_pipeline.lib.tasks.write_new_variants_table.load_gencode_gene_symbol_to_gene_id',
     )
-    @patch(
-        'v03_pipeline.lib.tasks.base.base_update_variant_annotations_table.UpdatedReferenceDatasetQueryTask',
-    )
     def test_sv_export_vcf(
         self,
-        mock_rd_query_task: Mock,
         mock_load_gencode: Mock,
     ) -> None:
         mock_load_gencode.return_value = GENE_ID_MAPPING
-        mock_rd_query_task.return_value = MockCompleteTask()
         worker = luigi.worker.Worker()
         update_variant_annotations_task = (
             UpdateVariantAnnotationsTableWithNewSamplesTask(