Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. 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.

How can I get variants only from homozygous region.

Hi,

I have DNASeq data of an offspring "C" which is cross of parents "A" and "B". Reference genomes of both "A" and "B" are available. I am interested in getting the variants, which should be in the "A" homozygous region. I am confused which reference to use and how to proceed. Your advice will be very valuable. Thank you.

Tagged:

Answers

  • tommycarstensentommycarstensen United KingdomMember ✭✭✭

    @bishwo, check out the --excludeIntervals (-XL) option available to all GATK Walkers for excluding variants, which are not homozygous:

    https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_engine_CommandLineGATK.php#--excludeIntervals
    

    You can use the flag multiple times. It is additive. You might also want to check --intervals and --interval_set_rule.

  • bishwobishwo Member

    @tommycarstensen, thanks for quick response. To which reference should I map C to, reference genome of A or that of B ?

  • tommycarstensentommycarstensen United KingdomMember ✭✭✭

    @bishwo These are not humans? You should use the same reference, which the reads were aligned to. You can check this in the bam header with samtools view -H your.bam. Check the SQ tags.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @bishwo‌

    Hi,

    Can you clarify on what exactly your end goal is?

    Thanks,
    Sheila

  • bishwobishwo Member

    @sheila I am interested in knowing "A" homozygous variants in the offspring "C". C is an offspring of "A" and "B".

    Firstly, I mapped "C" to the reference genome of "A" and called variants using HaplotypeCaller, but I did not get a single homozygous reference variants (0/0).

    Secondly, I mapped "C" to reference genome of "B" and I got quite many homozygous alternate (1/1) variants.

    I am wondering why I did not get any homozygous reference variants ("A" homozygous ) in the first run, but got homozygous alternate ("A" homozygous) variants in the second run ?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @bishwo‌

    Hi,

    If you used Haplotype Caller in regular mode, you will not get calls that are homozygous reference. Haplotype Caller will only output sites that have a variant allele, such as heterozygous or homozygous variant. So, all the sites that are not present in the output are homozygous reference.

    To get the sites that are homozygous reference, you can use HaplotypeCaller in GVCF mode, then use GenotypeGVCFs with -allSites.
    https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_haplotypecaller_HaplotypeCaller.php#--emitRefConfidence
    https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_GenotypeGVCFs.php#--includeNonVariantSites

    Thanks,
    Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    To add to Sheila's answer -- your analysis design is not appropriate for the question you want to answer. To use GATK, you should pick one individual as the reference and map all your samples to that reference. In a general case, you could pick either one of the parents or, better, another individual outside the trio that is considered representative of the population. In your case, one of the parents may be more appropriate than the other, but it depends what you mean by "A homozygous variants in C". Do you mean you want to know which sites in C are homozygous for an allele present in A but not in B? If so, the simplest way to do it is use B as reference; call variants in A and C, and select variants that share an ALT allele and are homozygous in C.

  • bishwobishwo Member

    Thank you @Sheila‌ and @Geraldine_VdAuwera‌ for you valuable comment.

    @Geraldine_VdAuwera‌ Yes, I want to know the sites in C that are homozygous for an allele present in A only. I have samples only from C. A and B are available as reference. However, only A has gene annotation available. Previously, I have mapped C samples to B reference and selected ALT alleles which are homozygous to C. Since I do not have gene annotation from B, I could not map the variants to genes. So, it seems the only way out is to map to A and get homozygous reference alleles using the method suggested by @Sheila‌.

  • bishwobishwo Member
    edited January 2015

    @Sheila @Geraldine_VdAuwera‌
    I have three C samples. I executed following command to to get sites that are homozygous reference.

    java -Xmx20g -jar $GATKPATH \ -T HaplotypeCaller \ -R $WHOLEGENOME \ -I $OutputDir/$sampleName-sorted-dupMarked-realigned-recalibrated.bam \ -o $OutputDir"-VCF"/$sampleName-raw-snp-HAPLOTYPECALLER.vcf \ --dbsnp $DBSNP \ -bamout $OutputDir/$sampleName-haplotype.bam \ --emitRefConfidence GVCF \ -variant_index_type LINEAR \ -variant_index_parameter 128000
    I used following command to run GenotypeGVCFs.

    java -Xmx20g -jar $GATKPATH \ -R $WHOLEGENOME \ -T GenotypeGVCFs \ --variant GATK-VCF/13_M-raw-snp-HAPLOTYPECALLER.vcf \ --variant GATK-VCF/16_M-raw-snp-HAPLOTYPECALLER.vcf \ --variant GATK-VCF/34_M-raw-snp-HAPLOTYPECALLER.vcf \ -allSites \ -o GATK-VCF/all.vcf
    However, I got following error when I tried to execute above GenotypeGVCFs.

    ##### ERROR Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
    ##### ERROR
    ##### ERROR MESSAGE: The list of input alleles must contain <NON_REF> as an allele but that is not the case at position 711; please use the Haplotype Caller with gVCF output to generate appropriate records
    ##### ERROR ------------------------------------------------------------------------------------------
    

    It was not recommend to post this error to forum, but I have run Haplotype Caller on gVCF mode.

    When I check 711th line I did not find NON_REF (as shown below). There are thousands of such lines without NON_REF as an allele.

    . . . Chr1 20485 . A <NON_REF> . . BLOCK_SIZE=10;END=20494 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:18:42:18:42 Chr1 20495 . A <NON_REF> . . BLOCK_SIZE=5;END=20499 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:19:36:18:33 Chr1 20500 . A <NON_REF> . . BLOCK_SIZE=9;END=20508 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:18:24:17:21 Chr1 20509 . AATAT A 193.73 . AC=1;AF=0.500;AN=2;BaseQRankSum=-0.555;ClippingRankSum=-0.660;DP=19;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=58.78;MQ0=0;MQRankSum=-0.312;QD=2.55;ReadPosRankSum=0.935 GT:AD:DP:GQ:PL 0/1:7,6:13:99:231,0,335 Chr1 20514 . A <NON_REF> . . BLOCK_SIZE=1;END=20514 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:19:0:19:0 Chr1 20515 . T <NON_REF> . . BLOCK_SIZE=31;END=20545 GT:DP:GQ:MIN_DP:MIN_GQ 0/0:16:45:15:42 . . .

    I could not figure out anything wrong in my script. I tried to search the meaning of NON_REF, but did not succeeded. What does it mean? Currently I am using GenomeAnalysisTK-3.3-0.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @bishwo Can you please confirm that what you are showing is the GVCF produced by HC without any modifications? And it was produced by GATK version 3.3-0?

  • bishwobishwo Member
    edited January 2015

    @Geraldine_VdAuwera‌ am sorry that I messed up with two versions in my pipeline. Now, I got GVCFs results. I found two different patters in the result files. Following is a part of output file that shows variants with no alternate base and variants with alternate base. How do i interpret the result ? How would I know if there is mutation in these sites ?

    Chr1    708     .       A       .       .       .       AN=6;DP=51;NCC=0        GT:DP   0/0:17  0/0:16  0/0:18
    Chr1    709     .       T       .       .       .       AN=6;DP=53;NCC=0        GT:DP   0/0:18  0/0:17  0/0:18
    Chr1    710     .       G       .       .       .       AN=6;DP=53;NCC=0        GT:DP   0/0:18  0/0:17  0/0:18
    Chr1    711     .       T       C       953.16  .       AC=4;AF=0.667;AN=6;BaseQRankSum=0.722;ClippingRankSum=0.844;DP=52;FS=0.000;GQ_MEAN=147.67;GQ_STDDEV=84.01;MLEAC=4;MLEAF=0.667;MQ=57.23;MQ0
    =0;MQRankSum=1.30;NCC=0;QD=18.33;ReadPosRankSum=0.722;SOR=0.643       GT:AD:DP:GQ:PL  0/1:8,10:18:99:246,0,189        0/1:9,8:17:99:203,0,231 1/1:0,17:17:51:536,51,0
    Chr1    712     .       T       .       .       .       AN=6;DP=53;NCC=0        GT:DP   0/0:18  0/0:17  0/0:18
    Chr1    713     .       T       .       .       .       AN=6;DP=54;NCC=0        GT:DP   0/0:18  0/0:18  0/0:18
    Chr1    714     .       G       .       .       .       AN=6;DP=52;NCC=0        GT:DP   0/0:18  0/0:16  0/0:18
    
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @bishwo‌

    Hi,

    If you have aligned your C samples to the reference of A, then the sites that are hom-ref (do not have an alternate allele) are A homozygous. The sites that have an alternate allele should have the allele from B at that position.

    Since you are looking for A homozygous, you can simply take the sites which have no alternate alleles.

    -Sheila

  • bishwobishwo Member

    @Sheila‌ Thank you very much for your reply. Yes, I have aligned C samples to the reference of A and got above result. Is there any way to find possible mutation on one of the alleles homozygous to A?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @bishwo‌

    Hi,

    Can you clarify what you mean by "possible mutation on one of the alleles homozygous to A"?

    This article may help you: http://gatkforums.broadinstitute.org/discussion/1268/how-should-i-interpret-vcf-files-produced-by-the-gatk

    -Sheila

  • bishwobishwo Member
    edited February 2015

    @Sheila Sorry for the confusion. Please ignore my previous comment.

    I aligned 3 C samples to the reference genome of B and selected homozygous ALT alleles (1/1). I filtered raw variants based on coverage cutoff 10 and quality cutoff 100. I got quite many variants from one of the samples i.e. 3_C. Following plot shows number of variants discovered among three samples. There is extremely huge number of variants in sample 3_C. All three samples are from same parents A and B. I am wondering if you have some insight regarding this variation.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @bishwo
    Hi,

    I cannot see the plot you posted. Are you saying there are too many variants in one sample compared to the other two? I am not sure I can help you answer that. You will have to do some quality control analysis on your sample data.

    -Sheila

Sign In or Register to comment.