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.

--snpmask in FastaAlternateReferenceMaker not inserting N?

hawchuanhawchuan WashingtonMember

Hi Sheila,
(PS I've posted this as a comment under another thread, but figured that the title of that other thread might be too obscure)
I read a few other posts regarding FastaAlternateReferenceMaker with snpmask. Here's how I constructed mine to mask sites (for one particular sample, just to test) that have low quality (DP below 8). I first emit all site (BP_RESOLUTION) under HC fand then use -allSites under GenotypeGVCFs. From this all-sites VCF, I then used SelectVariants to select sites that have DP < 8 under INFO. This presumably gives me a vcf-formated masking file for use in the final step. However, after I have run FastaAlternateReferenceMaker with this file and --snpmask, masked sites are still showing reference nucleotide, instead of N. Here's my last command:

java -jar GenomeAnalysisTK.jar \
-T FastaAlternateReferenceMaker \
-R arach.uce.mar17.cs.fa \
-o Al_12620.alt.ref.fasta \
-V Al_12620.allsites.vcf \
--snpmask mask.vcf \
--use_IUPAC_sample Index_Al_12620

What is it that I'm doing wrong here? I'm running version 3.3.
Many thanks!

HC

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @hawchuan
    Hi HC,

    Can you confirm that the sites not being called N are indeed SNPs in your mask.vcf and not no-called sites?

    Thanks,
    Sheil

  • hawchuanhawchuan WashingtonMember

    Hi Sheila,
    Thanks for the response. I followed a few of the discussions on this issue, including one started by bsarver (where u put in a request to call ./. as N rather than to skip over). I generated my all-sites vcf to be masked against using the following command:

    -T GenotypeGVCFs \
    -R myref \
    --variant my.gVCF \
    -allSites \
    -o Al_12620.allsites.vcf

    The result of this is that all sites with at least 1 read is called (ie 0/0. 0/1 or 1/1). Those sites with DP=0 is not called (ie ./.).
    I then made another vcf based on this one using SelectVariants to create a "mask" (see cmd below). My understanding is that, when applied to FastaAlternateReferenceMaker, this will force gatk to output N when a called site coincides with the a masked site.

    -R myref \
    -T SelectVariants \
    -V Al_12620.allsites.vcf \
    -select "DP < 8" \
    -o mask.vcf

    I checked against Al_12620.allsites.vcf and found that indeed called sites (ie not ./.) have been given reference base (instead of N) even when the sites are specified in the mask.

    For example, in Al_12620.allsites.vcf, the 2nd base of 1st locus is "called" as A (although DP =2). The first base is N because the ref base is an ambiguity.

    CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Index_Al_12620

    uce-7824_(modified)__cs 1 . N . . . AN=2;DP=2;NCC=0 GT:DP 0/0:2
    uce-7824_(modified)__cs 2 . A . . . AN=2;DP=2;NCC=0 GT:DP 0/0:2

    In the mask.vcf file, both sites are included by SelectVariant because both have DP below 8.

    CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Index_Al_12620

    uce-7824_(modified)__cs 1 . N . . . AN=2;DP=2;NCC=0 GT:DP 0/0:2
    uce-7824_(modified)__cs 2 . A . . . AN=2;DP=2;NCC=0 GT:DP 0/0:2

    Finally, the fasta ouput for the 1st locus is this

    1

    NA...
    when the 2nd base should be N instead of A (the ref) because it is specified by the mask.vcf

    I hope the above clarify somethings, not make it worse. Thanks!

  • SheilaSheila Broad InstituteMember, Broadie admin

    @hawchuan
    Hi HC,

    It looks like you are right. I'm not sure if the other users who have asked about this found the same problem, but I am going to add the issue to the original bug report (./. not being masked)

    Thanks,
    Sheila

  • hawchuanhawchuan WashingtonMember

    Hi Sheila,
    Thank you. That's a bummer. I was hoping it was just my mistake and it will be an easy fix. Will bug fixes normally come out with new releases (ie 3.5 being the next one) or put out as patches?
    HC

  • SheilaSheila Broad InstituteMember, Broadie admin

    @hawchuan
    Hi HC,

    I can't really predict when the fix will come out. It depends on when the developers get to the bug. I will see what I can do to get this moved up in priority, but there are some other bugs in front of yours that take precedence. If you would like, you can always submit a patch to us with the code fix, and we can merge it. That always takes less time.

    -Sheila

Sign In or Register to comment.