To celebrate the release of GATK 4.0, we are giving away free credits for running the GATK4 Best Practices pipelines in FireCloud, our secure online analysis portal. It’s first come first serve, so sign up now to claim your free credits worth $250. Sponsored by Google Cloud. Learn more at https://software.broadinstitute.org/firecloud/documentation/freecredits

GenotypeGVCFs sample specific genotype columns

Hi, I have generated 9 .g.vcf files from my .bam files, and want to do joint genotype calling with all these gvcf files using GenotypeGVCFs tool. I have the following command:

java -jar GenomeAnalysisTK.jar -T GenotypeGVCFs -R ref.fasta -V:VCF gvcflist.list -o output.vcf

In my output, I expected to find 9 sample columns for the genotype (ie 0/1:195,32:250:99:534,0,6705), however I can only find one sample column in the vcf file, and it seems that all 9 samples were pooled together. May I know if I have done something wrong? I would like to know individual genotype for each sample.

I have tried to use:

java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R ref.fasta -I sam1.bam -I sam2.bam -I sam3.bam ....

and also only obtain one column instead of multiple-sample column for the genotype.

Thank you.

Best Answer

Answers

  • scwong85scwong85 Member
    edited October 2015

    Hi Sheila,

    Thank you very much! You are absolutely right, I have the same read group for all my samples it seems :\ Is there anyway for me to change the read group in the .g.vcf file? Or do I need to use picard to replace my read group at the bam level and do the HaplotypeCaller again?

    Thank you.

    Edit: Yes I can just modify the "SM" value from the gvcf file and do GenotypeGVCFs on all the gvcf files in the list, and I will have multiple genotype columns. Thanks again Sheila for your help! :)

    Post edited by scwong85 on
  • thuthuynguyenthuthuynguyen ParisMember
    edited October 2015

    Hi Sheila,
    I am biologist and absolutely new to GATK and command lines. I have tried all Best Practices from GATK for my project in searching minor variants in HIV gen and I am stuck in the GenotypeVCFs step. I have several g.vcf files and create 1 vcf file as output by GenotypeVCFs command line but all samples were pooled together. This is the same problem with SCWONG85 as discussed here but I don't know how to fix it. I don't understand your response. How can I change the sample name in my read group or modify the "SM" value from the gvcf file? Can you tell me the command line to follow or an example for that. Sorry for the inconvenience. I really appreciate all of your advices for new beginner like me.
    Thanks so much
    Thuy NGUYEN

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @thuthuynguyen
    Hi Thuy Nguyen,

    There are two ways to do this.
    1) You can edit the bam header using Picard's AddOrReplaceReadGroups. http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups

    2) You can simply change the sample name in your GVCF. This is the easiest way. You can simply use a text editor to manually change the sample name to whatever you want.

    I hope this helps.

    -Sheila

  • nancySEEnancySEE malaysiaMember

    Firstly, i would like to thanks GATK team again for creating such a perfect feature to help to pool the samples together!
    My savior!!

    However, i do have a little issue with the combined result.

    I attempt to pool two samples (they are same different, but were sequenced with two different approaches: WGS and GBS methods) together from two gvcf files. However, my expected outcome is partially correct.

    #

    Example of a SNP loci in Sample1_WGS:
    LG1 81482 . T C, 149.77 . BaseQRankSum=-1.858;ClippingRankSum=0.454;DP=19;MLEAC=1,0;MLEAF=0.500,0.00;MQ=50.00;MQ0=0;MQRankSum=-0.372;ReadPosRankSum=0.289 GT:AD:DP:GQ:PL:SB 0/1:11,8,0:19:99:178,0,327,211,350,561:9,2,6,2

    Example of a SNP loci in Sample1_GBS:
    LG1 81482 . T C, 73.77 . BaseQRankSum=0.919;ClippingRankSum=-0.778;DP=14;MLEAC=1,0;MLEAF=0.500,0.00;MQ=50.00;MQ0=0;MQRankSum=-0.354;ReadPosRankSum=-0.071 GT:AD:DP:GQ:PL:SB 0/1:10,4,0:14:99:102,0,282,132,294,426:6,4,4,0

    Example of combined result obtained (Run options: -T GenotypeGVCFs; GATK version 3.3-0):
    LG1 81482 . T C 149.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.919;ClippingRankSum=0.454;DP=33;FS=0.000;GQ_MEAN=178.00;MLEAC=1;MLEAF=0.500;MQ=50.00;MQ0=0;MQRankSum=-3.540e-01;NCC=0;QD=7.88;ReadPosRankSum=0.289;SOR=0.399 GT:AD:DP:GQ:PL 0/1:11,8:19:99:178,0,327

    #

    The combined DP in INFO field seems correct, however looking at GT:AD:DP:GQ:PL, it report the result of first sample instead of reporting the combined result. Same goes to the quality field. Do you have any idea whether i had miss up anything? i had changed the sample name in the header of gvcf file.

    Looking forward to hear from you. =)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Glad you find the tools useful, @nancySEE!

    My guess is that the tool is not properly recognizing that there are two samples. Can you please post the headers of the two GVCF files?

  • nancySEEnancySEE malaysiaMember

    Many thanks Geraldine.
    Attached with the the headers of the two gvcf files. =)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Ah, see the last line of the header: both samples are still named Sample_17_BLN. You need to give them different names in those fields. Even if you change the names elsewhere, it is that place that matters.

  • nancySEEnancySEE malaysiaMember

    Hi Geradine, i would want pool to these two samples into the one column only.
    The reason behind is i do have another combined gvcf file that these two samples are in different columns (Eg: CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample_17_BLN_GBS Sample1_17_BLN_WGS)", thus i would want to pool these two same samples (different methods) as one (Eg: CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample_17_BLN), ultimately aim to compare the genotypes among three different approaches (GBS method only ; WGS method only ; combined of GBS and WGS methods)
    So i'm expected something like this 0/1:21,12:33:99:178,0,327 =)

    Issue · Github
    by Sheila

    Issue Number
    308
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • nancySEEnancySEE malaysiaMember

    For the combined of GBS and WGS approach, can i combine the two gvcf files straight (change the sample name to the same to do the genotypeGVCFs, so that the AD from two samples will be merged, and the PL & GT will be recalculated)? Is this a correct way to do this? Originally i' was thinking of merging the two BAM files into one, and redo the variant calling. Just wondering if there's any quickest alternative, for example combine the gvcf files (since i already have them at the first place).
    Cheers.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Oh I see -- that's not at all the intended use of the tools. If you're trying to see what the calls would look like if you pooled the data from these samples together, you should run them through HaplotypeCaller together with the same sample identifier (SM). Any other way of combining the calls after HaplotypeCaller will produce inferior results. I understand what you're trying to do but I'm afraid there is currently no tool in GATK that would apply the correct combination logic.

  • nancySEEnancySEE malaysiaMember

    Seems that there's no way i can go for shortcuts.
    Alright, i'll take your advice and re-perform HaplotypeCaller from the merged BAM files.
    Thanks again.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Actually I found out yesterday that we have an internal tool that does what you want, but it's currently not public. I'm checking if we can make it public and include it in the 3.5 release that is coming next week. No promises but I'll see what I can do.

  • nancySEEnancySEE malaysiaMember

    Such an exciting news for me!!!
    Looking forward to GATK new release with fingers crossed. =)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @nancySEE, sorry for getting back to you so late on this. After a bit of digging I found that the tool I mentioned is specifically designed for combining WGS + exome, and that there is still some development work that needs to happen before it can be applied to other use cases. I'm afraid we won't be able to release this in the upcoming 3.5 version, sorry to disappoint you! But I will let you know when it is ready.

  • nancySEEnancySEE malaysiaMember

    Hi Geraldine, I knew you had did your very best!
    Will keep an eye on the upcoming release if this feature integrated.
    Cheers and thanks for the effort!

  • I have the exact same problem as scwong85, but I checked the read groups in the BAM files, and they all have different names. In the mapping stage, I gave the reads names based on the sample name and a lane identifier. So what's wrong?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @olavur
    Hi,

    Sorry for the delay. Can you tell us if the GVCF files have different sample names in them? If so, which sample name is present in the final VCF?

    Thanks,
    Sheila

  • olavurolavur Member
    edited August 2017

    @Sheila thanks for the reply.

    I may have found my error, I'll try it and get back to you.

    Post edited by olavur on
  • @Sheila it turned out to be a user error, just a dumb mistake. Sorry about that.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @olavur
    Hi,

    Thanks for letting us know. Glad you solved it :smile:

    -Sheila

Sign In or Register to comment.