SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
calling variants using GATK kjaja Bioinformatics 0 03-19-2013 10:39 AM
SNP calling using GATK UnifiedGenotyper baika Bioinformatics 1 02-25-2013 01:41 PM
Problems using GATK for variant calling pooja_singh RNA Sequencing 1 07-07-2012 10:23 AM
GATK snp calling wanguan2000 Bioinformatics 0 11-24-2011 08:23 PM
SNP calling around junction sites chenjy Bioinformatics 2 07-29-2011 11:12 PM

Reply
 
Thread Tools
Old 05-02-2014, 02:03 AM   #1
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default Calling invariant sites with GATK

Hello,

I would like to obtain a vcf with variant AND INVARIANT sites using GATK.
I am not sure about the options to use in order to obtain invariant sites.
I tried a few but still I have only variants.

Anyone can help ?

Thanks a lot !

Muriel
MurielGB is offline   Reply With Quote
Old 05-02-2014, 06:12 AM   #2
vdauwera
Member
 
Location: Boston, MA

Join Date: Apr 2012
Posts: 42
Default

Hi Muriel,

What you want is to run the GATK's HaplotypeCaller in GVCF mode, with the arguments --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000 added to your command line. Have a look at the documentation here: http://www.broadinstitute.org/gatk/g...-discovery-ovw and here: http://www.broadinstitute.org/gatk/g...rticle?id=3893

These documents cover joint analysis of cohorts, but the basic principle is the same for single samples, you just skip the joint genotyping step.
vdauwera is offline   Reply With Quote
Old 05-02-2014, 08:28 AM   #3
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

Hi,
Thank you vdauwera for your answer.
The thing is that when you do that, you can't have a bam containing multiple samples.
Muriel
MurielGB is offline   Reply With Quote
Old 05-02-2014, 08:48 AM   #4
cariboudoug
Member
 
Location: Toronto

Join Date: Oct 2013
Posts: 17
Default

If you want to use GATK's UnifiedGenotyper you just need to add: --output EMIT_ALL_SITES
cariboudoug is offline   Reply With Quote
Old 05-02-2014, 08:53 AM   #5
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

Thank you cariboudoug !
Indeed, I tried this already but I was thinking that it's probably better to use HaplotypeCaller since this is the new version of UnifiedGenotyper...
MurielGB is offline   Reply With Quote
Old 05-02-2014, 09:26 AM   #6
vdauwera
Member
 
Location: Boston, MA

Join Date: Apr 2012
Posts: 42
Default

Quote:
Originally Posted by MurielGB View Post
Hi,
Thank you vdauwera for your answer.
The thing is that when you do that, you can't have a bam containing multiple samples.
Muriel
Oh, if you have multiple samples, you can just call each in turn using the GVCF mode, then you do the joint genotyping to get a single multi-sample VCF with the results.

If your problem is that you have all the samples in one bam, the simplest is to split it into per-sample bams, or you can pass the bam with a read filter to specify a single sample to run on in each run.
vdauwera is offline   Reply With Quote
Old 05-04-2014, 12:52 AM   #7
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

Hi !

I had a look on this new pipeline.

I launched variant calling on a bam file containing only one sample with the options mentioned before :
Quote:
java -Xmx300g -jar /GenomeAnalysisTK/3.0.0/GenomeAnalysisTK.jar \
-T HaplotypeCaller -R ref.fasta \
-I ${INPUT}.bam \
--genotyping_mode DISCOVERY \
-stand_emit_conf 0 \
-stand_call_conf 0 \
-o ${INPUT}\_VC2.vcf \
--emitRefConfidence GVCF \
--variant_index_type LINEAR \
--variant_index_parameter 128000 \
-nct 16
I still don't have all positions in the output vcd as you can see :
Quote:
KE340296.1 822 . C <NON_REF> . . END=823GTP:GQ:MIN_DP:PL 0/0:7:18:7:0,18,270
KE340296.1 824 . C <NON_REF> . . END=824GTP:GQ:MIN_DP:PL 0/0:7:0:7:0,0,219
KE340296.1 825 . C <NON_REF> . . END=827GTP:GQ:MIN_DP:PL 0/0:7:18:7:0,18,270
KE340296.1 828 . G <NON_REF> . . END=829GTP:GQ:MIN_DP:PL 0/0:7:0:7:0,0,199
KE340296.1 830 . T <NON_REF> . . END=854GTP:GQ:MIN_DP:PL 0/0:6:18:6:0,12,180
KE340296.1 855 . C <NON_REF> . . END=855GTP:GQ:MIN_DP:PL 0/0:7:1:7:0,1,233
KE340296.1 856 . T <NON_REF> . . END=890GTP:GQ:MIN_DP:PL 0/0:4:6:3:0,6,90
KE340296.1 891 . G <NON_REF> . . END=905GTP:GQ:MIN_DP:PL 0/0:4:0:0:0,0,0
KE340296.1 906 . T <NON_REF> . . END=944GTP:GQ:MIN_DP:PL 0/0:2:6:2:0,6,37
KE340296.1 945 . C <NON_REF> . . END=104GTP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0
Do you know why ?

Muriel
MurielGB is offline   Reply With Quote
Old 05-04-2014, 08:27 AM   #8
vdauwera
Member
 
Location: Boston, MA

Join Date: Apr 2012
Posts: 42
Default

Sure, if you want a call for every single position (instead of averaged blocks, which are used to compress the information in order to minimize file size) you just change -ERC GVCF to -ERC BP_RESOLUTION.
vdauwera is offline   Reply With Quote
Old 05-04-2014, 10:48 PM   #9
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

Oh great !!
I was wondering what was the difference between these two options !
Thanks a lot for your help vdauwera !
MurielGB is offline   Reply With Quote
Old 05-05-2014, 12:49 AM   #10
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

Well... I am wondering if I won't have another problem using -ERC BP_RESOLUTION.
Because after that I obtain a vcf file per sample containing both variant and invariant sites, I want to run the joint analysis to have a vcf combining all my samples.

In the documentation regarding the GenotypeGVCF option that perform the joint analysis (here : https://www.broadinstitute.org/gatk/...typeGVCFs.html), it is written that : GenotypeGVCFs merges gVCF records that were produced as part of the "single sample discovery" pipeline using the '-ERC GVCF' mode of the Haplotype Caller.

So I am wondering if this is gonna work since I use -ERC BP_RESOLUTION rather than GVCF ?

Thanks a lot !

Muriel
MurielGB is offline   Reply With Quote
Old 05-05-2014, 02:24 PM   #11
vdauwera
Member
 
Location: Boston, MA

Join Date: Apr 2012
Posts: 42
Default

That's completely fine actually, the BP_RESOLUTION gvcfs are also a valid input for the joint genotyping step. I'll update the documentation to be more accurate and inclusive. Feel free to post any further questions or comments about GATK tools on our support forum
vdauwera is offline   Reply With Quote
Old 05-05-2014, 10:01 PM   #12
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

So cooool !!!
Thanks a lot, I am very happy of this new tool, it's gonna allow population geneticists to do so much !
MurielGB is offline   Reply With Quote
Old 05-06-2014, 07:01 AM   #13
MurielGB
Member
 
Location: Montpellier, France

Join Date: Oct 2013
Posts: 51
Default

I have another problem, sorry vdauwera

I did variant calling on two samples individually using the option described and it works well, I have both variant and invariant sites.

However, the QUAL and INFO fields are missing in most cases :
Quote:
KE332545.1 64 . G <NON_REF> . . . GT:ADP:GQ:PL 0/0:17,0:17:51:0,51,660
KE332545.1 65 . G <NON_REF> . . . GT:ADP:GQ:PL 0/0:17,0:17:51:0,51,655
KE332545.1 66 . G T,<NON_REF> 1.14 . BaseQRankSum=-2.602;ClippingRankSum=-0.903;DP=24;MLEAC=1,0;MLEAF=0.500,0.00;MQ=38.57;MQ0=0;MQRankSum=1.115;ReadPosRankSum=0.372 GT:ADP:GQ:PL:SB 0/1:14,4,0:18:23:23,0,449,64,460,524:10,4,3,1
KE332545.1 67 . T <NON_REF> . . . GT:ADP:GQ:PL 0/0:18,0:18:53:0,53,583
KE332545.1 68 . T <NON_REF> . . . GT:ADP:GQ:PL 0/0:18,0:18:54:0,54,659
KE332545.1 69 . T <NON_REF> . . . GT:ADP:GQ:PL 0/0:18,0:18:54:0,54,699
It seems that these two fields are present only in case of an alternative allele existing right ?
Why ? How I can I use my filtering pipeline for the invariant sites !

Then I performed the joint analysis and this is the same.
This is a shame because I really hoped I will also have these informations for invariant sites...
Is there something that I missed ?

Thank you,

Muriel
MurielGB is offline   Reply With Quote
Old 10-17-2014, 11:13 PM   #14
rach16
Junior Member
 
Location: SF bay area

Join Date: Oct 2014
Posts: 2
Default

Hi Muriel

Were you able to find an answer to your question?

Thanks
rach16 is offline   Reply With Quote
Old 06-21-2017, 02:03 AM   #15
gshieh2
Junior Member
 
Location: Taiwan

Join Date: Jun 2017
Posts: 2
Default

We need the "QUAL" of invariant sites (homogeneous reference sites) from GATK. Can we just plug "0" into the formula -10 log_10 Prob (0 variant) since ALT is ". "?
gshieh2 is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:53 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO