GenotypeGVCFs can't handle somatic GVCFs from DRAGEN 3.8.4

[ldgauthier]

Bug Report

Affected tool(s) or class(es)

GenotypeGVCFs

Affected version(s)

• Latest public release version [version?]
• Latest master branch as of [date of test?]

Description

GenotypeGVCFs won't joint call DRAGEN mitochondrial data because of the DRAGEN somatic output format.

We should be able to use the DRAGEN SQ in place of Mutect2's TLOD (see line 279 in GenotypeGVCFsEngine)
Note that DRAGEN SQ is a Phred-scaled double

Steps to reproduce

DRAGEN somatic GVCF entries from version 3.8.4 look like:
chrM 1 . G <NON_REF> . weak_evidence END=1 GT:AD:DP:SQ:MIN_DP 0/0:112,1579:1691:0:1691

Run GenotypeGVCFs with -V to a file like that (reference GenotypeGVCFsIntegrationTest::testGenotypingForSomaticGVCFs() for more details)
Must include --input-is-somatic as of now

Expected behavior

The task should run to completion, calculating a site quality store using the DRAGEN SQ value

Actual behavior

Error from AlleleFrequencyCalculator about not having PLs or GQ.

[ldgauthier]

Replaces #7797

[MartaRusmini]

Hi
I'm trying to merge three VCF files from a trio but I have a problem with the chrM. I generated the gvcf files with DRAGEN (v. 4.2.4) germline calling.
My starting VCF file contains the SQ detail and it is like this:
chrM 263 rs2853515 A G <NON_REF> . PASS DP=1216;MQ=245.48;FractionInformativeReads=0.978;SoftClipRatio=0;DB GT:SQ:AD:AF:F1R2:F2R1:DP:SB:MB:PS 1|1:98.12 0:0 1173 0:1 0:0 368 0:0 805 0:1173:0 0 306 867:0 0 599 574:263

After merging, I obtain a file with all the information except the genotype:
chrM 263 rs2853515 A G,<NON_REF> . . DP=3248;FractionInformativeReads=0.995;MQ=247.66;SoftClipRatio=0.00 GT:AD:AF:DP:F1R2:F2R1:MB:PS:SB:SQ .|.:1,1198,0:0.9992,0.0000:1199:0,471,0:1,727,0:0,1,607,591:263:1,0,286,912:98.12,0.00 .|.:0,1173,0:1.0000,0.0000:1173:0,368,0:0,805,0:0,0,599,574:263:0,0,306,867:98.12,0.00 .|.:0,824,0:1.0000,0.0000:824:0,326,0:0,498,0:0,0,393,431:263:0,0,171,653:98.12,0.00

I used the following command:
gatk CombineGVCFs -O TEST.gvcf.gz --variant first.hard-filtered.gvcf.gz --variant /second.hard-filtered.gvcf.gz --variant third.hard-filtered.gvcf.gz -R /hg38/Homo_sapiens_assembly38.fasta --intervals chrM

After reading your post, I also tried
gatk CombineGVCFs -O somat.gvcf.gz --input-is-somatic --variant first.hard-filtered.gvcf.gz --variant /second.hard-filtered.gvcf.gz --variant third.hard-filtered.gvcf.gz -R /hg38/Homo_sapiens_assembly38.fasta --intervals chrM
but I obtained the same result.

Am I wrong?
To be sure, I also tried to use GenomicDBImport, but the list of variants was completely empty.

Can you help me please?
thanks
Marta

[ldgauthier]

Hi @MartaRusmini -- the typical CombineGVCFs behavior is to output no-call genotypes, like you see here, and then GenotypeGVCFs reassigns the genotypes using the supplied (or default) quality threshold. I'll take a look and see if it's a quick fix to implement the suggestion I had years ago

[ldgauthier]

Hi @MartaRusmini -- I may have a solution that works. Would you be willing to share the whole chrM GVCF for the three samples?

[MartaRusmini]

Hi @ldgauthier
Sorry for the delay and thank you for your help!
Here, the three chrM gvcf files

T-NGS_9008.chrM.gvcf.gz
T-NGS_9009.chrM.gvcf.gz
T-NGS_9010.chrM.gvcf.gz

Thank you very much

[ldgauthier]

Perfect timing, thanks @MartaRusmini ! I have a PR with some very simple test data (#9188) and running on this data will make me much more confident of the new functionality. We're hoping to get this merged and do a release this month.