Half-bad hetvars: one bad allele fails the genotype [VS-1615] [VS-1649] #9149

mcovarr · 2025-04-14T23:11:31Z

VETS VCFs

Look for differences between the expected and actual VCFs for VETS with a command like this:

diff <(grep -E -v '^#' actual/0000000000-quickit.vcf) <(grep -E -v '^#' expected/0000000000-quickit.vcf) | grep -A 3 high_CALIBRATION_SENSITIVITY_SNP

VETS SNP example: new behavior on the upper line, old behavior on the lower line. Allele 1 has a calibration sensitivity of 0.9977, greater than the 0.997 threshold. With the new behavior, all three samples fail. The first two samples have 0/1 genotypes and the third has 1/2. With the old behavior, the first two samples fail but the third is rescued by the 2 allele's low calibration sensitivity.

< chr20 26078548        .       A       G,T     .       .       AS_QUALapprox=0|1642|245;CALIBRATION_SENSITIVITY=0.9977,0.7783;QUALapprox=546;SCORE=-0.6604,-0.4697    GT:AD:FT:GQ:RGQ 0/1:15,13,0:high_CALIBRATION_SENSITIVITY_SNP:99:511     0/1:13,14,0:high_CALIBRATION_SENSITIVITY_SNP:99:546    1/2:0,13,5:high_CALIBRATION_SENSITIVITY_SNP:99:756
---
> chr20 26078548        .       A       G,T     .       .       AC=1,1;AF=0.500,0.500;AN=2;AS_QUALapprox=0|1642|245;CALIBRATION_SENSITIVITY=0.9977,0.7783;QUALapprox=546;SCORE=-0.6604,-0.4697 GT:AD:FT:GQ:RGQ 0/1:15,13,0:high_CALIBRATION_SENSITIVITY_SNP:99:511    0/1:13,14,0:high_CALIBRATION_SENSITIVITY_SNP:99:546     1/2:0,13,5:PASS:99:756

VETS INDEL example: Here allele 2 has a calibration sensitivity of 0.993, greater than the INDEL threshold of 0.99. In the upper line (new behavior), sample 3 with genotype 1/2 fails. On the lower line (old behavior), sample 3 is rescued by the relatively low calibration sensitivity of allele 1.

< chr20 1882863 .       CTTTTTTTTT      C,CTTTTTTTTTTT,CTTTTTTTTTT      .       .       AC=1,0,1;AF=0.250,0.00,0.250;AN=4;AS_QU
ALapprox=0|1034|59|146;CALIBRATION_SENSITIVITY=0.8956,0.993,0.941;QUALapprox=526;SCORE=-0.5033,-0.6215,-0.5311  GT:AD:FT:GQ:PGT
:PID:RGQ        0/1:7,14,0,0:PASS:99:1|0:1882818_A_T:567        0/3:7,0,0,9:PASS:88:.:.:146     1/2:0,12,2,0:high_CALIBRATION_S
ENSITIVITY_INDEL:59:.:.:526
---
> chr20 1882863 .       CTTTTTTTTT      C,CTTTTTTTTTTT,CTTTTTTTTTT      .       .       AC=2,1,1;AF=0.333,0.167,0.167;AN=6;AS_Q
UALapprox=0|1034|59|146;CALIBRATION_SENSITIVITY=0.8956,0.993,0.941;QUALapprox=526;SCORE=-0.5033,-0.6215,-0.5311 GT:AD:GQ:PGT:PI
D:RGQ   0/1:7,14,0,0:99:1|0:1882818_A_T:567     0/3:7,0,0,9:88:.:.:146  1/2:0,12,2,0:59:.:.:526
7825c7825

VQSR VCFs

VQSR SNP example: new behavior on the upper line, old behavior on the lower line. The low_VQSLOD_SNP threshold here is -2.6332. Allele 1 has a VQSLOD below the threshold, but in the upper line the 1/2 genotype of the second and third samples PASSes. In the lower line this genotype fails.

< chr20 30218456        .       C       T,A     .       .       AC=2,2;AF=0.500,0.500;AN=4;AS_QUALapprox=0|2519|1809;AS_VQSLOD=
-6.4503,0.1166;QUALapprox=1072  GT:AD:FT:GQ:PGT:PID:RGQ 0/1:30,10,0:low_VQSLOD_SNP:99:0|1:30218454_G_A:326      1/2:14,26,14:PA
SS:68:.:.:1072  1/2:0,30,36:PASS:99:.:.:2495
> chr20 30218456        .       C       T,A     .       .       AS_QUALapprox=0|2519|1809;AS_VQSLOD=-6.4503,0.1166;QUALapprox=1
072     GT:AD:FT:GQ:PGT:PID:RGQ 0/1:30,10,0:low_VQSLOD_SNP:99:0|1:30218454_G_A:326      1/2:14,26,14:low_VQSLOD_SNP:68:.:.:1072
        1/2:0,30,36:low_VQSLOD_SNP:99:.:.:2495

VQSR INDEL example: new behavior on the upper line, old behavior on the lower line. The low_VQSLOD_INDEL threshold here is -0.5365. Allele 1 has a VQSLOD below the threshold, but in the upper line the 4/1 genotype of the third samples PASSes. In the lower line this genotype fails.

< chr20 817544  .       GAAAAAA GAAAAAAA,G,GAAAAA,GA    .       .       AC=1,1,1,1;AF=0.250,0.250,0.250,0.250;AN=4;AS_QUALapprox=0|111|481|111|246;AS_VQSLOD=-1.4797,5.3448,3.8604,3.5201;QUALapprox=275       GT:AD:FT:GQ:RGQ 2/3:0,0,12,5,0:PASS:99:582    4/1:0,2,0,0,6:PASS:29:275        0/1:1,5,0,0,0:low_VQSLOD_INDEL:2:82
---
> chr20 817544  .       GAAAAAA GAAAAAAA,G,GAAAAA,GA    .       .       AC=0,1,1,0;AF=0.00,0.500,0.500,0.00;AN=2;AS_QUALapprox=0|111|481|111|246;AS_VQSLOD=-1.4797,5.3448,3.8604,3.5201;QUALapprox=275 GT:AD:FT:GQ:RGQ 2/3:0,0,12,5,0:PASS:99:582      4/1:0,2,0,0,6:low_VQSLOD_INDEL:29:275  0/1:1,5,0,0,0:low_VQSLOD_INDEL:2:82

VDS Creation

This simply uses the new(ish) GvsTieOutVDS to tie out new VETS VCFs with a newly created VDS. Link to run here.

VDS Merge & Rescore

The changes to the merge & rescore logic are exactly the same as those in the import logic. Here we perform a merge and rescore run with this updated 1/2 logic, and then crack open the output VDS to spot check.

First check the input "echo" VDS and see that a particular hetvar genotype with one bad allele is not filtered (FT flag is True despite the failing "2" allele):

>>> def filter(vd):
...   fvd = vd.filter_rows((vd.locus.contig == 'chr20') & (vd.locus.position == 1882863))
...   ffvd = fvd.filter_cols(fvd.s == 'ERS4367795')
...   return ffvd
...
>>> e = filter(echo_vds.variant_data)
>>> e.entries().show()
+---------------+-------------------------------------------------+----------+                                     (0 + 1) / 1]
| locus         | alleles                                         | filters  |
+---------------+-------------------------------------------------+----------+
| locus<GRCh38> | array<str>                                      | set<str> |
+---------------+-------------------------------------------------+----------+
| chr20:1882863 | ["CTTTTTTTTT","C","CTTTTTTTTTT","CTTTTTTTTTTT"] | {}       |
+---------------+-------------------------------------------------+----------+

+---------------------------------------------------------------------------------------------+--------------+-------+-------+
| as_vets                                                                                     | s            |    GQ |   RGQ |
+---------------------------------------------------------------------------------------------+--------------+-------+-------+
| dict<str, struct{model: str, calibration_sensitivity: float64}>                             | str          | int32 | int32 |
+---------------------------------------------------------------------------------------------+--------------+-------+-------+
| {"C":("INDEL",8.50e-01),"CTTTTTTTTTT":("INDEL",9.21e-01),"CTTTTTTTTTTT":("INDEL",9.95e-01)} | "ERS4367795" |    59 |   526 |
+---------------------------------------------------------------------------------------------+--------------+-------+-------+

+-------+------+--------------+--------------+------+
|    PS | LGT  | LAD          | LA           |   FT |
+-------+------+--------------+--------------+------+
| int64 | call | array<int32> | array<int32> | bool |
+-------+------+--------------+--------------+------+
|    NA | 1/2  | [0,12,2]     | [0,1,3]      | True |
+-------+------+--------------+--------------+------+

Now looking at the "output" VDS that has gone through merge and rescore, a genotype with this same failing allele has FT of False:

>>> o = filter(out_vds.variant_data)
>>> o.entries().show()
+---------------+-------------------------------------------------+----------+                                     (0 + 1) / 1]
| locus         | alleles                                         | filters  |
+---------------+-------------------------------------------------+----------+
| locus<GRCh38> | array<str>                                      | set<str> |
+---------------+-------------------------------------------------+----------+
| chr20:1882863 | ["CTTTTTTTTT","C","CTTTTTTTTTT","CTTTTTTTTTTT"] | {}       |
+---------------+-------------------------------------------------+----------+

+---------------------------------------------------------------------------------------------+--------------+-------+-------+
| as_vets                                                                                     | s            |    GQ |   RGQ |
+---------------------------------------------------------------------------------------------+--------------+-------+-------+
| dict<str, struct{model: str, calibration_sensitivity: float64}>                             | str          | int32 | int32 |
+---------------------------------------------------------------------------------------------+--------------+-------+-------+
| {"C":("INDEL",5.60e-01),"CTTTTTTTTTT":("INDEL",9.33e-01),"CTTTTTTTTTTT":("INDEL",9.92e-01)} | "ERS4367795" |    59 |   526 |
+---------------------------------------------------------------------------------------------+--------------+-------+-------+

+-------+------+--------------+--------------+-------+
|    PS | LGT  | LAD          | LA           |    FT |
+-------+------+--------------+--------------+-------+
| int64 | call | array<int32> | array<int32> |  bool |
+-------+------+--------------+--------------+-------+
|    NA | 1/2  | [0,12,2]     | [0,1,3]      | False |
+-------+------+--------------+--------------+-------+

…1615]

…r_vcfs

mcovarr · 2025-04-22T16:00:31Z

.dockstore.yml

@@ -173,7 +173,7 @@ workflows:
       branches:


All the changes in this file were just for testing purposes and don't need to be merged.

mcovarr · 2025-04-22T16:01:12Z

scripts/variantstore/scripts/vds_validation.py

@@ -17,34 +17,34 @@



this is all stowaway changes fixing some bonkers whitespace that was driving IntelliJ nuts

mcovarr · 2025-04-22T16:01:32Z

scripts/variantstore/wdl/test/GvsQuickstartIntegration.wdl

@@ -178,7 +178,7 @@ workflow GvsQuickstartIntegration {
                git_branch_or_tag = git_branch_or_tag,
                git_hash = GetToolVersions.git_hash,
                use_VETS = false,
-                extract_do_not_filter_override = true,
+                extract_do_not_filter_override = false,


will revert before merge

mcovarr · 2025-04-22T16:02:00Z

scripts/variantstore/wdl/test/MakeHetvarVDS.wdl

@@ -0,0 +1,105 @@
+version 1.0


don't think this needs to be merged, revert

gbggrant · 2025-04-23T17:43:30Z

scripts/variantstore/scripts/import_gvs.py

-            any_indel_ok=acc.any_indel_ok | (~allele_is_snp[called_idx - 1] & allele_OK[called_idx - 1]),
-        ), hl.struct(any_no=False, any_yes=False, any_snp=False, any_indel=False, any_snp_ok=False, any_indel_ok=False)))
+            all_snps_ok=acc.all_snps_ok & (~allele_is_snp[called_idx - 1] | allele_OK[called_idx - 1]),
+            all_indels_ok=acc.all_indels_ok & (allele_is_snp[called_idx - 1] | allele_OK[called_idx - 1]),


I am confused as to why this is not:

Suggested change

all_indels_ok=acc.all_indels_ok & (allele_is_snp[called_idx - 1] | allele_OK[called_idx - 1]),

all_snps_ok=acc.all_snps_ok & (allele_is_snp[called_idx - 1] & allele_OK[called_idx - 1]),

all_indels_ok=acc.all_indels_ok & (~allele_is_snp[called_idx - 1] & allele_OK[called_idx - 1]),

Maybe I'm just confounded by the logic

Taking the SNP line as an example, I interpreted this as "nothing wrong SNP-wise if it's not a SNP at all, OR (implying it is a SNP) if the allele is OK".

If the allele was bad and it was an INDEL, the following line would catch it.

Actually now that I look at this again I think I can/should clean it up a bit... I don't think we care about "any" SNPs or INDELs any more now that the default for the "ok" accumulators is True.

…r_vcfs

mcovarr added 16 commits April 14, 2025 18:39

Half-bad hetvars, VCF edition: one bad allele fails the genotype [VS-…

f75c299

…1615]

dockstore

1d30881

VQSR version

e21f0ac

VDS attempt

3ee277a

Merge remote-tracking branch 'origin/ah_var_store' into vs_1615_hetva…

756d363

…r_vcfs

dockstore

4e18503

make VDS from half-bad failed hetvar VCFs

0171302

update merge and rescore too

d5b88b6

turn on filtering for VQSR manual tieout

d428e13

Merge remote-tracking branch 'origin/ah_var_store' into vs_1615_hetva…

182019c

…r_vcfs

variants docker

815c1a3

dockstore

bda44b7

holy moly tabs

45b65b2

docker

ed09ccc

dockstore

ace768f

dockstore

852e96f

mcovarr marked this pull request as ready for review April 22, 2025 15:52

mcovarr changed the title ~~Half-bad hetvars, VCF edition: one bad allele fails the genotype [VS-1615]~~ Half-bad hetvars: one bad allele fails the genotype [VS-1615] [VS-1649] Apr 22, 2025

mcovarr commented Apr 22, 2025

View reviewed changes

remove stuff prior to review

7b47b30

gbggrant reviewed Apr 23, 2025

View reviewed changes

mcovarr added 3 commits April 23, 2025 14:14

Merge remote-tracking branch 'origin/ah_var_store' into vs_1615_hetva…

7c9e753

…r_vcfs

docker

eaf3b08

include Hail optimizations, update Docker

f919d24

gbggrant approved these changes Apr 24, 2025

View reviewed changes

mcovarr added 2 commits April 24, 2025 11:35

simplify more

1a21244

Docker

23d45a6

mcovarr added 28 commits April 25, 2025 13:16

another hash bump

a545738

Merge remote-tracking branch 'origin/ah_var_store' into vs_1615_hetva…

fdb1988

…r_vcfs

how did that get messed up

2311f27

spanning deletion . quality scores become NaN doubles

7dc7c85

fresh baked GATK Docker with . quality fixes

ab94c99

hash bump for all chr run

b3bb025

hash bump

218dae4

hash bump

41f3844

hash bump

27ac2bb

hash bump

8e52897

vat hash bump

0f21f6f

vat hash bump

6f831cb

vat hash bump

64a5bcc

vat hash bump

3a9a0ae

hash bump

6f606f0

hash bump

347eb2f

hash bump

6d92d22

Merge remote-tracking branch 'origin/ah_var_store' into vs_1615_hetva…

521a606

…r_vcfs

cleanup

5fbc434

hash bump

a8e5796

hash bump

33a67a8

hash bump

2a8a8b5

hash bump

e2bc338

cleanup, VAT from VDS

09a0d5f

hash bump

6c6fa2b

hash bump

9ade1a4

hash bump

c09d1a3

hash bump

7aa2ce1

mcovarr merged commit 8c93b77 into ah_var_store May 8, 2025
20 checks passed

mcovarr deleted the vs_1615_hetvar_vcfs branch May 8, 2025 22:06

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Half-bad hetvars: one bad allele fails the genotype [VS-1615] [VS-1649] #9149

Half-bad hetvars: one bad allele fails the genotype [VS-1615] [VS-1649] #9149

Uh oh!

mcovarr commented Apr 14, 2025 •

edited

Loading

Uh oh!

mcovarr Apr 22, 2025

Uh oh!

mcovarr Apr 22, 2025

Uh oh!

mcovarr Apr 22, 2025

Uh oh!

mcovarr Apr 22, 2025

Uh oh!

gbggrant Apr 23, 2025

Uh oh!

mcovarr Apr 23, 2025 •

edited

Loading

Uh oh!

mcovarr Apr 23, 2025

Uh oh!

Uh oh!

Uh oh!

	all_indels_ok=acc.all_indels_ok & (allele_is_snp[called_idx - 1] \| allele_OK[called_idx - 1]),
	all_snps_ok=acc.all_snps_ok & (allele_is_snp[called_idx - 1] & allele_OK[called_idx - 1]),
	all_indels_ok=acc.all_indels_ok & (~allele_is_snp[called_idx - 1] & allele_OK[called_idx - 1]),

Half-bad hetvars: one bad allele fails the genotype [VS-1615] [VS-1649] #9149

Half-bad hetvars: one bad allele fails the genotype [VS-1615] [VS-1649] #9149

Uh oh!

Conversation

mcovarr commented Apr 14, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

VETS VCFs

VQSR VCFs

VDS Creation

VDS Merge & Rescore

Uh oh!

mcovarr Apr 22, 2025

Choose a reason for hiding this comment

Uh oh!

mcovarr Apr 22, 2025

Choose a reason for hiding this comment

Uh oh!

mcovarr Apr 22, 2025

Choose a reason for hiding this comment

Uh oh!

mcovarr Apr 22, 2025

Choose a reason for hiding this comment

Uh oh!

gbggrant Apr 23, 2025

Choose a reason for hiding this comment

Uh oh!

mcovarr Apr 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mcovarr Apr 23, 2025

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

mcovarr commented Apr 14, 2025 •

edited

Loading

mcovarr Apr 23, 2025 •

edited

Loading