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

How to mark specific individual genotype as No Call

Hi,

I have been trying to find a way to mark specific individual genotypes as No Call.

I know that in VariantFiltration it is possible to add the option --setFilteredGTToNocall in order to mark filtered genotypes as no call. However, in my case, there is no available filter corresponding to my criteria. Let me explain. I have some diploid-haploid paired related samples, i.e. I expect the haploid individual to share one of the two alleles from its diploid related. Therefore, for positions where there are discordant genotypes, I would like to mark individual genotypes as as no call in order to not include potentials errors in my data. I couldn't find a way to apply this rational to the pipeline of VariantFiltration or SelectVariants...

Eventually, I manage to manipulate the .vcf file by myself to mark any discordant genotypes as . or ./. However, when later I want to select the positions for which there are maximum of 0.2 ratio of no call (using the option --maxNOCALLfraction of SelectVariants from the nightly build version), there was an error message "there aren't enough comumns for line...". I suppose that it's because I did manipulate the vcf file by myself and there could have been errors (although I'm quite confident with my script).

Therefore, I think it's better to not manipulate the vcf file and try to do things using the pipeline....but how?

Issue · Github
by Sheila

Issue Number
898
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
vdauwera

Best Answer

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Hi @nkobmoo,

    You're right that manipulating the VCF file directly often causes formatting errors that break the file. The VCF is really hard to work with that way. So let's try to solve this in a more robust way. However I'm not sure I understand exactly what you describe. Would you be able to post a few examples of positions where you know that the conditions you are looking for are fulfilled, and which genotype(s) should be set as filtered?

  • nkobmoonkobmoo ParisMember
    edited September 2016

    Hi @Geraldine_VdAuwera ,

    For examples, there are samples for which I expect to have the same genotypes at certain sites. Whenever this is not the case, I want to mark the genotypes at those sites as NoCall.

    For example:
    s
    caffold1 159775 . C G 1356.24 PASS AC=1;AF=0.333;AN=3;DP=82;ExcessHet=3.0103;FS=0.000;MQ=60.00;QD=29.75;SOR=2.376 GT:AD:DP:GQ:PL 0/0:51,0:51:99:0,108,2120 1:0,31:31:99:1395,0 [followed by other samples.....]

    As you can see, the genotypes of the first and the second sample are discordant. Therefore I want to mark these genotypes as NoCall (./. and . respectively).

    I couldn't find a way to do this using GATK pipeline. Finally, I managed to modifly my vcf file using a home-made script.

    Cheers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @nkobmoo
    Hi,

    There really is no easy way to set those sites to all no-calls. Have a look at this article under "Using VariantContext directly". You will need to check for all samples being hom-ref or het or hom-var with VariantFiltration. Then, you can use SelectVariants to set those sites to no-calls.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @nkobmoo
    Hi again,

    I realized I misunderstood this question as I was answering your other question.

    I have two questions:
    1) Are all your samples paired (one diploid and one haploid paired together)?
    2) Do you need the sites to be no-calls or can they be completely excluded from the output?

    Thanks,
    Sheila

    P.S. What kind of data are you working with? It sounds interesting :smile:

  • nkobmoonkobmoo ParisMember

    @Sheila

    To answer to your questions:

    1. Some samples are paired (one diploid sample paired with one haploid sample). Some others are not (only haploid samples without paired diploid or vice-versa).

    2. I would need the sites to be no-calls only for the discordant genotypes because some other paired samples may have concordant genotypes and I don't want to exclude these informations.

    For information, I'm working with fungi. I extract spores from a reproductive structure for dna extraction. Usually, I could do single-spore extractions, making those samples haploid. However, sometimes single-spore extraction is not possible. Therefore, I had to proceed with multiple-spores extractions, taken as proxies for diploids. In some case cases, I could do from the same reproductive structure both single-spore extractions and multiple-spore extractions, making them as paired samples......All this sound a bit complicated but it's the only way I can optimize my sampling.

    Issue · Github
    by Sheila

    Issue Number
    1282
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
Sign In or Register to comment.