Merge pull request #1050 from broadinstitute/dev

Dev

Author: jklugherz

v03_pipeline/lib/annotations/sv.py

@@ -180,6 +180,7 @@ def gnomad_svs(
     gnomad_svs_ht: hl.Table,
     **_: Any,
 ) -> hl.Expression:
+    gnomad_svs_ht = gnomad_svs_ht.drop('locus', 'alleles')
     return gnomad_svs_ht.annotate(
         ID=gnomad_svs_ht.KEY,
     )[ht['info.GNOMAD_V4.1_TRUTH_VID']]

v03_pipeline/lib/reference_datasets/gnomad_svs.py

@@ -14,5 +14,4 @@ def get_ht(path: str, reference_genome: ReferenceGenome) -> hl.Table:
         N_HET=ht.info.N_HET,
         N_HOMREF=ht.info.N_HOMREF,
     )
-    ht = ht.key_by('KEY')
-    return ht.drop('locus', 'alleles')
+    return ht.key_by('KEY')

v03_pipeline/lib/reference_datasets/gnomad_svs_test.py

@@ -0,0 +1,93 @@
+import unittest
+from unittest.mock import patch
+
+import hail as hl
+
+from v03_pipeline.lib.model import ReferenceGenome
+from v03_pipeline.lib.reference_datasets.reference_dataset import ReferenceDataset
+
+TEST_GNOMAD_SVS_RAW_HT = (
+    'v03_pipeline/var/test/reference_datasets/raw/gnomad_svs_from_vcf.ht'
+)
+
+
+class GnomadSVsTest(unittest.TestCase):
+    @patch('v03_pipeline.lib.reference_datasets.gnomad_svs.vcf_to_ht')
+    def test_gnomad_svs(self, mock_vcf_to_ht):
+        mock_vcf_to_ht.return_value = hl.read_table(TEST_GNOMAD_SVS_RAW_HT)
+        ht = ReferenceDataset.gnomad_svs.get_ht(ReferenceGenome.GRCh38)
+        self.assertEqual(
+            ht.collect(),
+            [
+                hl.Struct(
+                    KEY='gnomAD-SV_v3_BND_chr1_1a45f73a',
+                    locus=hl.Locus(
+                        contig='chr1',
+                        position=10434,
+                        reference_genome=ReferenceGenome.GRCh38,
+                    ),
+                    alleles=['N', '<BND>'],
+                    AF=0.11413399875164032,
+                    AC=8474,
+                    AN=74246,
+                    N_HET=8426,
+                    N_HOMREF=28673,
+                ),
+                hl.Struct(
+                    KEY='gnomAD-SV_v3_BND_chr1_3fa36917',
+                    locus=hl.Locus(
+                        contig='chr1',
+                        position=10440,
+                        reference_genome=ReferenceGenome.GRCh38,
+                    ),
+                    alleles=['N', '<BND>'],
+                    AF=0.004201000090688467,
+                    AC=466,
+                    AN=110936,
+                    N_HET=466,
+                    N_HOMREF=55002,
+                ),
+                hl.Struct(
+                    KEY='gnomAD-SV_v3_BND_chr1_7bbf34b5',
+                    locus=hl.Locus(
+                        contig='chr1',
+                        position=10464,
+                        reference_genome=ReferenceGenome.GRCh38,
+                    ),
+                    alleles=['N', '<BND>'],
+                    AF=0.03698499873280525,
+                    AC=3119,
+                    AN=84332,
+                    N_HET=3115,
+                    N_HOMREF=39049,
+                ),
+                hl.Struct(
+                    KEY='gnomAD-SV_v3_BND_chr1_933a2971',
+                    locus=hl.Locus(
+                        contig='chr1',
+                        position=10450,
+                        reference_genome=ReferenceGenome.GRCh38,
+                    ),
+                    alleles=['N', '<BND>'],
+                    AF=0.3238990008831024,
+                    AC=21766,
+                    AN=67200,
+                    N_HET=21616,
+                    N_HOMREF=11909,
+                ),
+                hl.Struct(
+                    KEY='gnomAD-SV_v3_DUP_chr1_01c2781c',
+                    locus=hl.Locus(
+                        contig='chr1',
+                        position=10000,
+                        reference_genome=ReferenceGenome.GRCh38,
+                    ),
+                    alleles=['N', '<DUP>'],
+                    AF=0.0019970000721514225,
+                    AC=139,
+                    AN=69594,
+                    N_HET=139,
+                    N_HOMREF=34658,
+                ),
+            ],
+        )

v03_pipeline/lib/reference_datasets/reference_dataset.py

@@ -104,6 +104,12 @@ def version(self, reference_genome: ReferenceGenome) -> str:
             )
         return version
 
+    def dataset_types(
+        self,
+        reference_genome: ReferenceGenome,
+    ) -> frozenset[DatasetType]:
+        return CONFIG[self][reference_genome][DATASET_TYPES]
+
     @property
     def enums(self) -> dict | None:
         return CONFIG[self].get(ENUMS)
@@ -143,11 +149,9 @@ def get_ht(
         if enum_selects:
             ht = ht.transmute(**enum_selects)
         ht = filter_contigs(ht, reference_genome)
-        # Reference Datasets are DatasetType agnostic, but these
-        # methods (in theory) support SV/GCNV.  SNV_INDEL
-        # is passed as a proxy for non-SV/GCNV.
-        validate_allele_type(ht, DatasetType.SNV_INDEL)
-        validate_no_duplicate_variants(ht, reference_genome, DatasetType.SNV_INDEL)
+        for dataset_type in self.dataset_types(reference_genome):
+            validate_allele_type(ht, dataset_type)
+            validate_no_duplicate_variants(ht, reference_genome, dataset_type)
         # NB: we do not filter with "filter" here
         # ReferenceDatasets are DatasetType agnostic and that
         # filter is only used at annotation time.

v03_pipeline/migrations/annotations/0004_add_gnomad_svs.py

@@ -0,0 +1,52 @@
+import hail as hl
+
+from v03_pipeline.lib.annotations import sv
+from v03_pipeline.lib.migration.base_migration import BaseMigration
+from v03_pipeline.lib.model import DatasetType, ReferenceGenome
+from v03_pipeline.lib.reference_datasets.reference_dataset import ReferenceDataset
+
+# This vcf was generated with the gatk command:
+#
+# gatk SVConcordance --verbosity DEBUG --evaluation /var/seqr/phase4.seqr.gnomad_v4_tmp.vcf.gz
+# --truth  /var/seqr/gnomad.v4.1.sv.sites.modified.vcf.bgz
+# --sequence-dictionary gs://gcp-public-data--broad-references/hg38/v0/Homo_sapiens_assembly38.dict
+#
+# Followed by:
+# bcftools annotate --rename-annots /var/seqr/remap /var/seqr/phase4.seqr.gnomad_v4_tmp.vcf.gz | bgzip > /var/seqr/phase4.seqr.gnomad_v4.vcf.gz
+#
+# where remap contains "INFO/TRUTH_VID GNOMAD_V4.1_TRUTH_VID"
+PHASE_4_CALLSET_WITH_GNOMAD_V4 = 'gs://seqr-loading-temp/phase4.seqr.gnomad_v4.vcf.gz'
+
+
+class AddGnomadSVs(BaseMigration):
+    reference_genome_dataset_types: frozenset[
+        tuple[ReferenceGenome, DatasetType]
+    ] = frozenset(
+        ((ReferenceGenome.GRCh38, DatasetType.SV),),
+    )
+
+    @staticmethod
+    def migrate(ht: hl.Table, **_) -> hl.Table:
+        mapping_ht = (
+            hl.import_vcf(
+                PHASE_4_CALLSET_WITH_GNOMAD_V4,
+                reference_genome=ReferenceGenome.GRCh38.value,
+                force_bgz=True,
+            )
+            .key_rows_by('rsid')
+            .rows()
+        )
+        ht = ht.annotate(
+            **{
+                'info.GNOMAD_V4.1_TRUTH_VID': mapping_ht[ht.key].info[
+                    'GNOMAD_V4.1_TRUTH_VID'
+                ],
+            },
+        )
+        gnomad_svs_ht = ReferenceDataset.gnomad_svs.get_ht(ReferenceGenome.GRCh38)
+        ht = ht.annotate(gnomad_svs=sv.gnomad_svs(ht, gnomad_svs_ht))
+        ht = ht.drop('info.GNOMAD_V4.1_TRUTH_VID')
+        return ht.annotate_globals(
+            versions=ht.globals.versions.annotate(gnomad_svs='1.0'),
+            enums=ht.globals.enums.annotate(gnomad_svs=hl.Struct()),
+        )

v03_pipeline/var/test/reference_datasets/GRCh38/gnomad_svs/1.0.ht/.README.txt.crc

@@ -1,3 +1,3 @@
 This folder comprises a Hail (www.hail.is) native Table or MatrixTable.
   Written with version 0.2.133-4c60fddb171a
-  Created at 2025/02/16 18:40:38
+  Created at 2025/03/05 12:27:53

v03_pipeline/var/test/reference_datasets/GRCh38/gnomad_svs/1.0.ht/.metadata.json.gz.crc

@@ -0,0 +1,3 @@
+This folder comprises a Hail (www.hail.is) native Table or MatrixTable.
+  Written with version 0.2.133-4c60fddb171a
+  Created at 2025/03/04 14:19:46