Question about QUAL in DeepTrio output

[shadrinams]

Hello,

I am running DeepTrio for a dataset with known true-positive SNPs and indels. I followed guidelines for DeepTrio but had to change --config DeepVariantWGS to DeepVariant_unfiltered at the glnexus_cli step as a default QUAL threshold of 10 removed a lot of my TP calls.
I have compared distributions of QUAL score in the TP subset and all calls found by DeepTrio. Please see an attached histogram of all DeepTrio calls vs TP calls.

Could you please tell me if it is expected that QUAL of TP calls is between 0 and30, while there are calls with QUAL of up to 100? If not, what I could do wrong?
Thank you!

Best regards,
Maria.

Setup

• Operating system: Linux
• DeepVariant version: deepvariant_deeptrio-1.1.0.sif
• Installation method (Docker, built from source, etc.): singularity/3.6.4
• Type of data: (sequencing instrument, reference genome, anything special that is unlike the case studies?) HiSeq X Ten, hg38, WGS

[AndrewCarroll]

Hi @shadrinams

From your description, it looks like you are running the DeepTrio merge component in the way that I would recommend.

Can I ask a few questions:

1. Are you taking the values from the QUAL field of the multisample glnexus VCF
2. How have you determined the TP sites? Are these Genome in a Bottle, or do they come from some other source.

Your QUAL values for all calls match expectations for DeepTrio. However, I would not expect that "true" variants would have a distribution which departs from all calls. I can think of one reason that this might be the case:

Are these true variants de novos? DeepTrio's quality distribution for de novo variants is very different from its general quality distribution. This occurs because DeepTrio has learned that de novo events are quite rare, and so requires a higher standard of evidence to make a call which is a de novo. In these cases, DeepTrio is not extremely confident in the call, which results in a lower quality value.

The other things it might be good to look at is the genotype quality (GQ) field. This is the most direct measure of DeepTrio's confidence in a call, and maps directly from the probability for the called class. The QUAL value of a multisample VCF comes from GLnexus, and is a bit less direct measure of call confidence (still, I don't expect the distributions to be off in the manner that you see).

I don't think I can immediately answer this puzzle, but hopefully with a little more information we can figure it out.

Thanks,
Andrew

[shadrinams]

Hello Andrew,

Thank you for your quick reply!

1. Are you taking the values from the QUAL field of the multisample glnexus VCF

Yes, exactly!

2. How have you determined the TP sites? Are these Genome in a Bottle, or do they come from some other source.

These are PCGC data, TPs were determined by combination of methods and manually curated. We expect an accuracy of the
found TPs to be > 95% (based on PCR for a similar dataset), although we might still miss some TP calls.

Are these true variants de novos? DeepTrio's quality distribution for de novo variants is very different from its general quality distribution. This occurs because DeepTrio has learned that de novo events are quite rare, and so requires a higher standard of evidence to make a call which is a de novo. In these cases, DeepTrio is not extremely confident in the call, which results in a lower quality value.

I am sorry, I did not mention it. Yes, we are looking for denovos in trios. We are comparing efficiency of a few methods to create a pipeline for a big dataset. I thought we might use the QUAL score from DeepTrio to filter calls found by GATK4 pipeline.
If we use GQ fields for further filtering what values do you recommend for parents and proband?

Now, I use the following filters to retrieve denovo calls from the multisample glnexus VCF:

• Heterozygous ratio of proband = 0.2-0.8
• Homozygous ratio of parents <= 0.1
• ALT allele depth of proband >= 7
• Genotype quality of proband >= 60
• Read depth >= 7
• Allele count = 1
• Some regional filters were applied to remove noisy regions
• Common variants were removed based on 1000genome and gnomad population frequencies
  Also, I had to split multiallelic calls and recalculate genotypes based on AD fields as I had a lot of ./. and 0/1 for Homozygous reference calls in parents.

As results, I obtained 909 SNPs and 1,236 indels for my 10 test trios. My list of TPs contains 698 SNPs and 61 indels. So, I still have a lot of false-positives calls. Is there a way to filter my variants from DeepTrio further?

Thank you!

Best regards,
Maria.

[AndrewCarroll]

Hi @shadriniams

Thank you for this additional information, it's very helpful.

Yes, for de novo variants, we do expect DeepTrio to be less confident in the variant position, so your observation makes sense.

The best field to filter on further is the genotype quality (likely of the proband). We haven't had large orthogonal datasets of de novo variants yet to give you a good recommendation for the threshold. Could you try to rank-order the calls by proband GQ and see if this separates the true and false calls reasonably well.

Finally, are you looking at sites that are 0/1 in child an 0/0 in both parents, or are you also looking at homozygous Mendelian violations. If it is the former, I don't think there's anything else to consider. However, if it is the latter, I might have some additional recommendations.

Thanks,
Andrew

[shadrinams]

Could you try to rank-order the calls by proband GQ and see if this separates the true and false calls reasonably well.

Please see the attached plot. It seems to tend to the same behavior. Looking at my plots I would think to filter denovo calls by the low part of QUAL or proband GQ scales ( QUAL/GQ < 30-40 ). Counterintuitive but seems pretty solid based on the plots.

Finally, are you looking at sites that are 0/1 in child an 0/0 in both parents, or are you also looking at homozygous Mendelian violations. If it is the former, I don't think there's anything else to consider. However, if it is the latter, I might have some additional recommendations.

I am looking for the former but we have plans to try Mendelian violations too. I would appreciate for your input for it!

Thank you!

Best regards,
Maria.

[AndrewCarroll]

Hi @shadrinams

Thank you for the plot. This is expected the the found de novo calls are lower in confidence (because DeepTrio has learned that de novo events are rare). Given that a call is a de novo (0/1-0/0-0/0), the higher GQ values will still indicate higher confidence, so more confident de novo calls should be more likely to be true.

For calls which are 1/1-0/0-1/1 (or permutations of this), the parent and child models of DeepTrio do not coordinate, so they aren't optimizing for consistency. There is a property of some regions which look like potential segmental duplications where a call that appears to a human as a 1/1/ or 0/1 is actually some kind of CNV. DeepTrio has learned some parts which are predictive of this (generally a variant-dense haplotype and a reference haplotype with higher depth). There may be cases where the child model gives a call a 1/1 and the parent gives a 0/0 when the site itself is similar, but the context is different. If you are looking for homozygous Mendelian violations, you may want to filter regions of high variant density, as apparent calls will come from this phenomenon.

In addition, if you are looking at the sex chromosomes, it would be good to separately call the non-PAR chrX for male samples. For a male sample, running chrX providing only the mother sample provides best results. (chrY should already only have paternal reads outside of the PAR).

Thank you,
Andrew

[shadrinams]

Hello Andrew,

Thank you for the plot. This is expected the the found de novo calls are lower in confidence (because DeepTrio has learned that de novo events are rare). Given that a call is a de novo (0/1-0/0-0/0), the higher GQ values will still indicate higher confidence, so more confident de novo calls should be more likely to be true.

I have tried to verify the idea and plotted normalized histograms for all DeepTrio calls, filtered DeNovo-like (0/1-0/0-0/0) calls and true DeNovo calls (verified in IGV). Please see an attached picture. True DeNovo calls are still in the middle-low part of proband GQ range for DeNovo-like calls. I have tried to look at several DeNovo-like variants with a higher proband GQ, which were not in my list of true DeNovos, but all of them were false-positives.

I have prepared a few examples of DeepTrio variants in IGV (squished and expanded views) with corresponding output in DeepTrio VCF to show you discrepancies between DeepTrio VCFs and BAMs. I think it could be a part of the GQ issue in DeNovos.
The order of FORMAT fields in multisample VCF is proband, mother, father.
The order of samples in IGV is father, mother, proband (from top to bottom).

1) True Denovo, QUAL=3, proband GQ=5. It looks good except very low proband GQ.

chr5 | 92696737 | chr5_92696737_C_T | C | T | 3 | . | AF=0.166667;AQ=3 | GT:DP:AD:GQ:PL:RNC | 0/1:32:17,15:5:3,0,32:.. | 0/0:27:27,0:50:0,108,1079:.. | 0/0:21:21,0:50:0,105,1049:..

**2) Filtered Denovo-like, QUAL=46, proband GQ=13. Mulitallelic, inherited; when VCF is normalized, it passes DeNovo filter as father's genotype misses ALT allele information **

chr5 | 24093912 | chr5_24093912_AAT_A;chr5_24093912_A_AATAT | AAT | A,AATATAT | 46 | . | AF=0.333333,0.166667;AQ=46,15 | GT:DP:AD:GQ:PL:RNC | 1/2:30:5,12,13:13:44,15,55,15,0,53:.. | 0/1:31:16,15,0:46:46,0,70,990,990,990:.. | ./.:30:27,1,0:18:0,18,45,990,990,990:II

chr5 | 24093912 | chr5_24093912_AAT_A;chr5_24093912_A_AATAT | AAT | AATATAT | 46 | . | AF=0.166667;AQ=15 | GT:DP:AD:GQ:PL:RNC | 0/1:30:5,13:13:44,15,53:.. | 0/0:31:16,0:46:46,990,990:.. | ./.:30:27,0:18:0,990,990:II

**3) Filtered Denovo-like, QUAL=27, proband GQ=28. Inherited; it passes DeNovo filter as father's genotype misses ALT allele information **

chr7 | 54624683 | chr7_54624683_A_AATC | A | AATC | 27 | . | AF=0.166667;AQ=27 | GT:DP:AD:GQ:PL:RNC | 0/1:39:22,16:28:27,0,48:.. | 0/0:40:40,0:50:0,120,1199:.. | 0/0:28:28,0:50:0,90,899:..

What do you think could be a reason for DeepTrio calling to miss some information from read allignments?
Thank you.

Best regards,
Maria.

[AndrewCarroll]

Hi @shadrinams

Your observations about some apparent de novos originating from undercalling of a parent is interesting. I think that we'll look into whether some harmonization or filtering after calls could make the resulting calls more uniform.

With respect to your examples. It is often difficult to be definitive about why DeepVariant does not make certain calls. I cannot give you a reason, but for some cases I can share observations. In all cases, I will list the call followed by observation.

chr5 | 92696737 | chr5_92696737_C_T | C | T | 3 | . | AF=0.166667;AQ=3 | GT:DP:AD:GQ:PL:RNC | 0/1:32:17,15:5:3,0,32:.. | 0/0:27:27,0:50:0,108,1079:.. | 0/0:21:21,0:50:0,105,1049:..

This looks clean. DeepTrio's GQ is low probably because it is a clear de novo and it has learned such events are rare.

chr5 | 24093912 | chr5_24093912_AAT_A;chr5_24093912_A_AATAT | AAT | A,AATATAT | 46 | . | AF=0.333333,0.166667;AQ=46,15 | GT:DP:AD:GQ:PL:RNC | 1/2:30:5,12,13:13:44,15,55,15,0,53:.. | 0/1:31:16,15,0:46:46,0,70,990,990,990:.. | ./.:30:27,1,0:18:0,18,45,990,990,990:II

One thing I note - it looks to me like there are 3 alleles represented in the reads for the top parent: 1) there is an insertion event in-phase with a downstream HET SNP. 2) There is a reference allele in-phase with REF at that later position. 3) There is evidence for a T SNP that is also in-phase with the downstream HET variant. For the reads that are HET T, it could be interesting to see if they overlap any other variants that would suggest that they come from a copy number variant elsewhere in the genome. It may be the case that DeepTrio does not call a variant in the parent because some of the variant reads may be coming from elsewhere.

chr7 | 54624683 | chr7_54624683_A_AATC | A | AATC | 27 | . | AF=0.166667;AQ=27 | GT:DP:AD:GQ:PL:RNC | 0/1:39:22,16:28:27,0,48:.. | 0/0:40:40,0:50:0,120,1199:.. | 0/0:28:28,0:50:0,90,899:..

This is interesting, since the evidence reported by DeepTrio doesn't match the view in the BAM. DeepTrio is saying that both parents don't have evidence for variant reads at this position. It is possible the after re-alignment (which DeepTrio does) there isn't evidence for a variant. However, the child is still called and has similar evidence. Is it possible for you to share snippets of this BAM file (e.g. a ~1000 bp window around this variant)? If so, I might like to take a look to see why DeepTrio reports no evidence in the reads.

Thanks,
Andrew

[shadrinams]

Hello Andrew,

Thank you for your reply.

I have prepared snippets of the BAM files for the trio with variant chr7:54624683. It would be interesting to see what you find!

One thing I note - it looks to me like there are 3 alleles represented in the reads for the top parent: 1) there is an insertion event in-phase with a downstream HET SNP. 2) There is a reference allele in-phase with REF at that later position. 3) There is evidence for a T SNP that is also in-phase with the downstream HET variant. For the reads that are HET T, it could be interesting to see if they overlap any other variants that would suggest that they come from a copy number variant elsewhere in the genome. It may be the case that DeepTrio does not call a variant in the parent because some of the variant reads may be coming from elsewhere.

It is an interesting thought! But even if there is an explanation for missing HET T, the AATATAT insertion allele is also missed in the father genotype. And for some reason it is shown as ./.

Best regards,
Maria.

chr7_54624683.tar.gz

[AndrewCarroll]

Hi @shadrinams

My apologies for the delay in reply. It took me awhile to get to this issue. I have some thoughts about how to improve consistency in the trio calling that may help some of these cases. But it may take a bit of time to explore those.

For the chr7_54624683 case, I suspect there is some interaction with the variant call and the normalization done by GLnexus. When I look at the VCF for the Proband and Father, I see:

Proband

chr7    54624686        .       A       ATC     61.6    PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:61:39:17,21:0.538462:61,0,7

Father

chr7    54624686        .       A       ATC     54.8    PASS    .       GT:GQ:DP:AD:VAF:PL      0/1:55:33:18,15:0.454545:54,0,68

In my GLnexus VCF, I also see the following line, which has a call in both the child and the father (the variant position is a bit different, but this does seem to reflect the reads I see in the sample).

chr7    54624686        chr7_54624686_A_ATC     A       ATC     61      .       AF=0.333333;AQ=61       GT:DP:AD:GQ:PL:RNC      0/1:33:18,15:55:54,0,68:..      0/0:32:32,0:50:0,180,1799:..    0/1:39:17,21:61:61,0,76:..

Would you mind pasting the command that you used to merge the gVCFs from this sample? When you jointly genotyped, were there also other samples at the time of joint genotyping, or was it only this sample?

For reference, this is the command that I used on the gVCFs of the files to generate this gVCF:

sudo apt-get -y install bcftools
sudo apt-get -y install tabix
sudo docker pull quay.io/mlin/glnexus:v1.2.7

sudo docker run \
     -v "${PWD}/outputs":"/outputs" \
     quay.io/mlin/glnexus:v1.2.7 \
     /usr/local/bin/glnexus_cli \
     --config DeepVariantWGS \
     /outputs/Child_chr7_54624683.g.vcf.gz \
     /outputs/Father_chr7_54624683.g.vcf.gz \
     /outputs/Mother_chr7_54624683.g.vcf.gz \
       | bcftools view - | bgzip -c > outputs/GitHub440.cohort.vcf.gz

[shadrinams]

Hello Andrew,

Thank you for the information! It is interesting!

chr7:54624686 A-ATC and chr7:54624683 A-AATC are different variants. I see both in my output.

But still my output for chr7:54624686 A-ATC is not the same. It was not called in my proband. The one in my father looks similar but QUAL and GQ are quite different. I think it comes from the DeepTrio calling step. How can our outputs be so different?
I use deepvariant_deeptrio-1.1.0.sif.

For any case, I show all three 54624683/54624686 variants here:
My proband:
chr7 54624683 . A AATC 27.8 PASS . GT:GQ:DP:AD:VAF:PL 0/1:28:39:22,16:0.410256:27,0,48
chr7 54624686 . A C 27.8 PASS . GT:GQ:DP:AD:VAF:PL 0/1:28:38:16,22:0.578947:27,0,53
chr7 54624686 . A ATC - no call

My father:
chr7 54624683 . A AATC - no call
chr7 54624686 . A C - no call
chr7 54624686 . A ATC 26.2 PASS . GT:GQ:DP:AD:VAF:PL 0/1:26:33:18,15:0.454545:26,0,44

My GLnexus VCF:
chr7 54624683 chr7_54624683_A_AATC A AATC 27 . AF=0.166667;AQ=27 GT:DP:AD:GQ:PL:RNC 0/1:39:22,16:28:27,0,48:.. 0/0:40:40,0:50:0,120,1199:.. 0/0:28:28,0:50:0,90,899:..
chr7 54624686 chr7_54624686_A_C;chr7_54624686_A_ATC A C 27 . AF=0.166667;AQ=27 GT:DP:AD:GQ:PL:RNC 0/1:38:16,22:28:27,0,53:.. 0/0:40:40,0:50:0,120,1199:.. 0/0:33:18,0:26:26,990,990:..
chr7 54624686 chr7_54624686_A_C;chr7_54624686_A_ATC A ATC 27 . AF=0.166667;AQ=26 GT:DP:AD:GQ:PL:RNC 0/0:38:16,0:28:27,990,990:.. 0/0:40:40,0:50:0,120,1199:.. 0/1:33:18,15:26:26,0,44:..

I use the following command for glnexus. I had to switch to --config DeepVariant_unfiltered as --config DeepVariantWGS filtered out most of my DeNovo calls as they tend to have QUAL of <20.

module load glnexus/1.2.7
glnexus_cli --config DeepVariant_unfiltered --threads $(nproc) $FATHER.g.vcf.gz $MOTHER.g.vcf.gz $CHILD.g.vcf.gz | bcftools view - | bgzip -c > ${FAMILY}.deeptrio.vcf.gz

Best regards,
Maria.

[rabiafidan]

Hi, I wondered if there are any updates on this issue. Were you able to separate the true positives from false positives?

[AndrewCarroll]

Hi @rabiafidan

We have found the following:

1. The use of GLnexus preset DeepVariant_unfiltered is preferred for retaining True de novo calls, and we have updated the documentation for this.
2. We also observe a reduction in 0/1 child, 0/0 parent calls when post-processing the final VCF to set a parent to ./. when that parent is 0/0, the child is 0/1 or 1/1, and that parent has either less than 8 reads covering the variant position, or and allele fraction of > 0.15.

The second filter seems to help reduces cases where there is not enough confidence to clearly call a de novo. Does this filtering strategy seem like it might further help refine your calls? We are considering whether to recommend postprocessing of this nature via a script in the future.

Thank you,
Andrew