bugfix: gcnv frequency annotations (#1016)

* Ensure gcnv frequency stats are updated on all variants present in callset

* Support annotating from the callset

* fix bad test

* add test

Author: bpblanken

v03_pipeline/lib/annotations/fields_test.py

@@ -199,7 +199,7 @@ def test_get_lookup_table_fields(
             list(
                 get_fields(
                     ht,
-                    DatasetType.SNV_INDEL.lookup_table_annotation_fns,
+                    DatasetType.SNV_INDEL.variant_frequency_annotation_fns,
                     lookup_ht=lookup_ht,
                     dataset_type=DatasetType.SNV_INDEL,
                     reference_genome=ReferenceGenome.GRCh38,

v03_pipeline/lib/annotations/gcnv.py

@@ -57,11 +57,11 @@ def end_locus(
     return hl.locus(ht.chr, ht.end, reference_genome.value)
 
 
-def gt_stats(ht: hl.Table, **_: Any) -> hl.Expression:
+def gt_stats(ht: hl.Table, callset_ht: hl.Table, **_: Any) -> hl.Expression:
     return hl.struct(
-        AF=hl.float32(ht.sf),
-        AC=ht.sc,
-        AN=hl.int32(ht.sc / ht.sf),
+        AF=hl.float32(callset_ht[ht.variant_id].sf),
+        AC=callset_ht[ht.variant_id].sc,
+        AN=hl.int32(callset_ht[ht.variant_id].sc / callset_ht[ht.variant_id].sf),
         Hom=hl.missing(hl.tint32),
         Het=hl.missing(hl.tint32),
     )

v03_pipeline/lib/model/dataset_type.py

@@ -277,7 +277,6 @@ def formatting_annotation_fns(
             ],
             DatasetType.GCNV: [
                 gcnv.end_locus,
-                gcnv.gt_stats,
                 gcnv.num_exon,
                 gcnv.sorted_gene_consequences,
                 gcnv.start_locus,
@@ -354,15 +353,18 @@ def genotype_entry_annotation_fns(self) -> list[Callable[..., hl.Expression]]:
         }[self]
 
     @property
-    def lookup_table_annotation_fns(self) -> list[Callable[..., hl.Expression]]:
+    def variant_frequency_annotation_fns(self) -> list[Callable[..., hl.Expression]]:
         return {
             DatasetType.SNV_INDEL: [
                 snv_indel.gt_stats,
             ],
             DatasetType.MITO: [
                 mito.gt_stats,
             ],
-        }[self]
+            DatasetType.GCNV: [
+                gcnv.gt_stats,
+            ],
+        }.get(self, [])
 
     @property
     def should_send_to_allele_registry(self):

v03_pipeline/lib/tasks/update_variant_annotations_table_with_deleted_families.py

@@ -45,7 +45,7 @@ def update_table(self, ht: hl.Table) -> hl.Table:
         ht = ht.annotate(
             **get_fields(
                 ht,
-                self.dataset_type.lookup_table_annotation_fns,
+                self.dataset_type.variant_frequency_annotation_fns,
                 lookup_ht=lookup_ht,
                 **self.param_kwargs,
             ),

v03_pipeline/lib/tasks/update_variant_annotations_table_with_deleted_project.py

@@ -43,7 +43,7 @@ def update_table(self, ht: hl.Table) -> hl.Table:
         ht = ht.annotate(
             **get_fields(
                 ht,
-                self.dataset_type.lookup_table_annotation_fns,
+                self.dataset_type.variant_frequency_annotation_fns,
                 lookup_ht=lookup_ht,
                 **self.param_kwargs,
             ),

v03_pipeline/lib/tasks/update_variant_annotations_table_with_new_samples.py

@@ -66,44 +66,42 @@ def update_table(self, ht: hl.Table) -> hl.Table:
 
         # Union with the new variants table and annotate with the lookup table.
         ht = ht.union(new_variants_ht, unify=True)
-        if self.dataset_type.has_lookup_table:
+        if self.dataset_type.variant_frequency_annotation_fns:
             callset_ht = get_callset_ht(
                 self.reference_genome,
                 self.dataset_type,
                 self.callset_path,
                 self.project_guids,
             )
+            lookup_ht = None
+            if self.dataset_type.has_lookup_table:
+                lookup_ht = hl.read_table(
+                    lookup_table_path(
+                        self.reference_genome,
+                        self.dataset_type,
+                    ),
+                )
+                # Variants may have fallen out of the callset and
+                # have been removed from the lookup table during modification.
+                # Ensure we don't proceed with those variants.
+                ht = ht.semi_join(lookup_ht)
+
             # new_variants_ht consists of variants present in the new callset, fully annotated,
             # but NOT present in the existing annotations table.
             # callset_variants_ht consists of variants present in the new callset, fully annotated,
             # and either present or not present in the existing annotations table.
             callset_variants_ht = ht.semi_join(callset_ht)
-            ht = ht.anti_join(callset_ht)
-            lookup_ht = hl.read_table(
-                lookup_table_path(
-                    self.reference_genome,
-                    self.dataset_type,
-                ),
-            )
+            non_callset_variants_ht = ht.anti_join(callset_ht)
             callset_variants_ht = callset_variants_ht.annotate(
                 **get_fields(
                     callset_variants_ht,
-                    self.dataset_type.lookup_table_annotation_fns,
-                    lookup_ht=hl.read_table(
-                        lookup_table_path(
-                            self.reference_genome,
-                            self.dataset_type,
-                        ),
-                    ),
+                    self.dataset_type.variant_frequency_annotation_fns,
+                    lookup_ht=lookup_ht,
+                    callset_ht=callset_ht,
                     **self.param_kwargs,
                 ),
             )
-            ht = ht.union(callset_variants_ht, unify=True)
-
-            # Variants may have fallen out of the callset and
-            # have been removed from the lookup table during modification.
-            # Ensure we don't proceed with those variants.
-            ht = ht.semi_join(lookup_ht)
+            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)

v03_pipeline/lib/tasks/update_variant_annotations_table_with_new_samples_test.py

@@ -48,10 +48,12 @@
 TEST_SNV_INDEL_VCF = 'v03_pipeline/var/test/callsets/1kg_30variants.vcf'
 TEST_SV_VCF = 'v03_pipeline/var/test/callsets/sv_1.vcf'
 TEST_GCNV_BED_FILE = 'v03_pipeline/var/test/callsets/gcnv_1.tsv'
+TEST_GCNV_BED_FILE_2 = 'v03_pipeline/var/test/callsets/gcnv_2.tsv'
 TEST_REMAP = 'v03_pipeline/var/test/remaps/test_remap_1.tsv'
 TEST_PEDIGREE_3 = 'v03_pipeline/var/test/pedigrees/test_pedigree_3.tsv'
 TEST_PEDIGREE_4 = 'v03_pipeline/var/test/pedigrees/test_pedigree_4.tsv'
 TEST_PEDIGREE_5 = 'v03_pipeline/var/test/pedigrees/test_pedigree_5.tsv'
+TEST_PEDIGREE_8 = 'v03_pipeline/var/test/pedigrees/test_pedigree_8.tsv'
 
 GENE_ID_MAPPING = {
     'OR4F5': 'ENSG00000186092',
@@ -1478,7 +1480,7 @@ def test_sv_update_vat(
             ],
         )
 
-    def test_gcnv_update_vat(
+    def test_gcnv_update_vat_multiple(
         self,
     ) -> None:
         worker = luigi.worker.Worker()
@@ -1614,3 +1616,58 @@ def test_gcnv_update_vat(
                 ),
             ],
         )
+        update_variant_annotations_task = (
+            UpdateVariantAnnotationsTableWithNewSamplesTask(
+                reference_genome=ReferenceGenome.GRCh38,
+                dataset_type=DatasetType.GCNV,
+                sample_type=SampleType.WES,
+                callset_path=TEST_GCNV_BED_FILE_2,
+                project_guids=['R0115_test_project2'],
+                project_remap_paths=['not_a_real_file'],
+                project_pedigree_paths=[TEST_PEDIGREE_8],
+                skip_validation=True,
+                run_id='second_run_id',
+            )
+        )
+        worker.add(update_variant_annotations_task)
+        worker.run()
+        ht = hl.read_table(update_variant_annotations_task.output().path)
+        self.assertEqual(ht.count(), 3)
+        self.assertCountEqual(
+            ht.select('gt_stats').collect(),
+            [
+                hl.Struct(
+                    # Existing variant, gt_stats are preserved
+                    variant_id='suffix_16456_DEL',
+                    gt_stats=hl.Struct(
+                        AF=hl.eval(hl.float32(4.401408e-05)),
+                        AC=1,
+                        AN=22720,
+                        Hom=None,
+                        Het=None,
+                    ),
+                ),
+                hl.Struct(
+                    # Exiting variant, gt_stats are updated
+                    variant_id='suffix_16457_DEL',
+                    gt_stats=hl.Struct(
+                        AF=8.111110946629196e-05,
+                        AC=3,
+                        AN=36986,
+                        Hom=None,
+                        Het=None,
+                    ),
+                ),
+                hl.Struct(
+                    # New variant.
+                    variant_id='suffix_99999_DEL',
+                    gt_stats=hl.Struct(
+                        AF=8.802817319519818e-05,
+                        AC=2,
+                        AN=22719,
+                        Hom=None,
+                        Het=None,
+                    ),
+                ),
+            ],
+        )

v03_pipeline/var/test/callsets/gcnv_2.tsv

@@ -0,0 +1,4 @@
+chr	start	end	name	sample	sample_fix	svtype	GT	CN	NP	QA	QS	QSE	QSS	ploidy	strand	variant_name	ID	rmsstd	defragmented	sf	sc	lt100_raw_calls	lt10_highQS_rare_calls	PASS_SAMPLE	PASS_FREQ	PASS_QS	HIGH_QUALITY	genes_any_overlap	genes_any_overlap_exonsPerGene	genes_any_overlap_totalExons	genes_strict_overlap	genes_strict_overlap_exonsPerGene	genes_strict_overlap_totalExons	genes_CG	genes_LOF	genes_any_overlap_Ensemble_ID	genes_LOF_Ensemble_ID	genes_CG_Ensemble_ID	var_source	callset_ovl	identical_ovl	partial_0_5_ovl	any_ovl	no_ovl	strvctvre_score	is_latest	sample_cram_basename
+chr1	100006937	100007881	COHORT_1_Y_cnv_16326	RP-2037_RGP_999_1_v1_Exome_GCP	RGP_999_1	DEL	1	1	1	4	4	4	4	2	-	suffix_16456	115884	0	FALSE	4.401408e-05	1	TRUE	TRUE	TRUE	TRUE	FALSE	FALSE	AC118553.2,SLC35A3	1,1	2	None	0	0	NA	AC118553.2,SLC35A3	ENSG00000283761.1,ENSG00000117620.15	ENSG00000283761.1,ENSG00000117620.15		round2	round1	cluster_6_CASE_cnv_74577	cluster_6_CASE_cnv_74577	cluster_6_CASE_cnv_74577	FALSE	0.583	FALSE	RGP_999_1
+chr1	100017586	100023212	CASE_10_X_cnv_38690	C1992_RGP_999_1_v1_Exome_C1992	RGP_999_1	DEL	1	0	2	20	30	30	20	2	-	suffix_16457	115887	0	FALSE	8.111111E-05	3	FALSE	FALSE	FALSE	TRUE	FALSE	FALSE	AC118553.2,SLC35A3	1,2	2	None	0	0	NA	AC118553.2,SLC35A3	ENSG00000283761.1,ENSG22222222222.15	ENSG00000283761.1,ENSG00000117620.15		round2	round1	cluster_22_COHORT_cnv_56835	cluster_22_COHORT_cnv_56835	cluster_22_COHORT_cnv_56835	FALSE	0.507	FALSE	RGP_999_1
+chr1	100017585	100023213	CASE_10_X_cnv_38690	C1992_RGP_999_1_v1_Exome_GCP	RGP_999_1	DEL	1	0	2	20	30	30	20	2	-	suffix_99999	115887	0	FALSE	8.802817e-05	2	FALSE	FALSE	FALSE	TRUE	FALSE	FALSE	AC118553.2,SLC35A3	1,2	3	None	0	0	NA	AC118553.2,SLC35A3	ENSG00000283761.1,ENSG00000117620.15	ENSG00000283761.1,ENSG00000117620.15		round2	round1	cluster_22_COHORT_cnv_56835	cluster_22_COHORT_cnv_56835	cluster_22_COHORT_cnv_56835	FALSE	0.507	FALSE	RGP_999_1

v03_pipeline/var/test/pedigrees/test_pedigree_8.tsv

@@ -0,0 +1,2 @@
+Project_GUID	Family_GUID	Family_ID	Individual_ID	Paternal_ID	Maternal_ID	Sex
+R0115_test_project2	family_3_1	family_3	RGP_999_1			F