Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

About VariantFiltration

monicafabbromonicafabbro ArgentinaMember

Hi Team!
I have selected those variants that are covered in my 7 samples (diploid organism) with AN=14 and then I applied the VariantFiltration walker in order to filter those that have less than 6X of coverage.

Command line:

Keep variants with AN=14 (7 samples / diploid organism = 7*2=14)

java -jar GenomeAnalysisTK.jar -T SelectVariants -R sequence.fasta -select "AN==14" --variant cohort.vcf -o cohort.AN.vcf

Mark variants having DP<6

java -jar GenomeAnalysisTK.jar -T VariantFiltration -R sequence.fasta --genotypeFilterExpression "DP<6" --genotypeFilterName LowDP --variant cohort.AN.vcf -o cohort.AN.DP.vcf

Remove LowDP marked variants

grep -vw "LowDP" cohort.AN.DP.vcf > cohort.filtered.vcf

But when I check this I see that in some variants my lowDP filter is not working, for example:
GL622784 1350823 . C T 8109.01 PASS AC=13;AF=0.929;AN=14;BaseQRankSum=-1.369e+00;ClippingRankSum=1.79;DP=176;FS=1.661;MLEAC=13;MLEAF=0.929;MQ=58.94;MQRankSum=2.78;QD=33.65;ReadPosRankSum=0.678;SOR=0.307 GT:AD:DP:GQ:PL
1/1:0,48:48:99:1|1:1350822_T_C:2247,156,0 1/1:0,22:22:72:1|1:1350797_C_T:1069,72,0 1/1:0,0:.:6:1|1:1350822_T_C:49,6,0
1/1:0,17:17:51:1|1:1350797_C_T:765,51,0 1/1:0,14:14:42:1|1:1350822_T_C:630,42,0 0/1:6,7:13:99:0|1:1350822_T_C:276,0,231
1/1:0,62:62:99:1|1:1350776_G_C:3106,211,0

Am I missing something here? I thought this should work.
Any hint on this will be very helpful

Thanks a lot!

MF

Answers

  • pdexheimerpdexheimer Member ✭✭✭✭

    Are all of your false negatives cases where the depth is completely missing? Note that even though it holds for DP, in general missing is not the same as zero - so VariantFiltration has to treat it separately. I think that it reacts by not considering that genotype for filtration, effectively leaving it unfiltered.

    You would need to specially test for DP being missing, but I'm afraid I don't know the best way to do that

  • monicafabbromonicafabbro ArgentinaMember

    Hi pedxheimer!
    Yes I have 169 cases where depth is missing like the example in my comment above.
    Why I am having a "." instead of a "0" in the DP field in the first place? What does it means? How I could treat that in order to filter them?
    My main goal here is to keep those variants with coverage in all my samples and to have at least 6X of coverage in all of them.
    Thanks in advance!

    MF

  • SheilaSheila Broad InstituteMember, Broadie admin

    @monicafabbro
    Hi MF,

    Unfortunately, Variant Filtration does not have an option to filter out . in DP. However, you are right the bigger question is why is there a . in DP when the AD is 0,0. You can try filtering on AD values. But, I think the example you posted above may be an indel issue. Are these calls from version 3.3 or lower?
    Can you please post the bam file and bamout file for the sample with this record: 1/1:0,0:.:6:1|1:1350822_T_C:49,6,0 ?

    Thanks,
    Sheila

  • monicafabbromonicafabbro ArgentinaMember

    Hi Sheila,
    Thanks for your reply.
    The calls are from the latest version of GATK.
    Do you need the entire bam file of this sample? Because it's 2,7Gb and BTW it's confidential data so if you could tell me where I can send it to you in a private way I would appreciate it.
    Thanks in advance!

    MF

  • monicafabbromonicafabbro ArgentinaMember

    Hi Sheila,
    Yes it stands for my name =).
    Here is the screenshot.
    Thanks a lot!!
    MF

  • SheilaSheila Broad InstituteMember, Broadie admin

    @monicafabbro
    Hi Monica,

    Are these the bamout files you posted? Is the sample with the 1/1:0,0:.:6:1|1:1350822_T_C:49,6,0 call at the bottom? Does this variant get called if you use -dontUseSoftClippedBases?

    Thanks,
    Sheila

  • monicafabbromonicafabbro ArgentinaMember

    Hi Sheila,
    Yes they are and yes the sample with the 1/1:0,0:.:6:1|1:1350822_T_C:49,6,0 call is at the bottom.
    I ran the Haplotype caller with the -dontUseSoftClippedBases parameter and still get that call:
    GL622784 1350823 . C T, 22.58 . DP=0;MLEAC=2,0;MLEAF=1.00,0.00;MQ=0.00 GT:GQ:PGT:PID:PL:SB 1/1:6:0|1:1350822_T_C:49,6,0,49,6,49:0,0,0,0
    Although, it has been called with a very low quality...

    MF

  • SheilaSheila Broad InstituteMember, Broadie admin

    @monicafabbro
    Hi Monica,

    Okay. It looks like there should be an insertion and/or a deletion in the region that would get rid of a lot of those SNPs. Can you confirm you used the GATK Best Practices pipeline? Did you do Indel Realignment and use Haplotype Caller?

    If so, I would like to have a look at this myself. Would you be able to upload a bug report? Instructions are here: http://gatkforums.broadinstitute.org/discussion/1894/how-do-i-submit-a-detailed-bug-report

    -Sheila

  • monicafabbromonicafabbro ArgentinaMember

    Hi Sheila,
    I am working with a non-standard reference so I did not run Indel realignment and Base Recalibration. I did run Haplotype caller with the ERC mode.
    How many lines of my BAM file would be suitable for you?

    Thanks in advance!
    MF

  • monicafabbromonicafabbro ArgentinaMember

    Hi Shelia,
    I have uploaded the bug report as you asked me with this file name: 20150607_gatk_depth_error.tar.gz
    Thanks

    MF

  • SheilaSheila Broad InstituteMember, Broadie admin

    @monicafabbro
    Hi Monica,

    I think this issue only arises when you run GenotypeGVCFs on many GVCFs. I cannot replicate the issue with only one sample. With one sample, the site of interest is not called. Can you upload one or two more samples that cause the issue to arise?

    Thanks,
    Sheila

  • monicafabbromonicafabbro ArgentinaMember

    Hi Sheila,
    I couldn't upload the file via sftp with the user and password GATK provide (I got an error message) but I uploaded the file using the user "anonymous".
    The file is in /incoming and it is called 20150729_gatk_depth_error.tar.gz.
    The file contains 7 samples including the one I uploaded before. That sample is called Sample07 now.
    Thank in advance!
    Monica

  • SheilaSheila Broad InstituteMember, Broadie admin

    @monicafabbro
    Hi Monica,

    I could not find your file. It seems anonymous does not work for uploading. The FTP should be back up and running, so you can try uploading again. Let me know if you still have issues.

    Thanks,
    Sheila

  • monicafabbromonicafabbro ArgentinaMember

    Hi Sheila,
    I have uploaded the file.
    Thanks!

    Mónica

  • SheilaSheila Broad InstituteMember, Broadie admin

    @monicafabbro
    Hi Monica,

    Because that site has a very shallow read depth, adding -minPruning 1 and -minDanglingBranchLength 1 makes the call have a proper DP. I think the low depth was throwing off the DP annotation for the site. I will make a note of this for future use.

    -Sheila

Sign In or Register to comment.