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 does VariantEval's CountVariants option determine MNPs?

Hi,

I am trying to understand what CountVariants is calling as an MNP in my dataset. When I run CountVariants, I get ~3K MNPs called across my samples. However, when I run SelectVariants with -selectType MNP (and -selectType MIXED), or when I run VariantAnnotator with -A VariantType, nothing seems to be identified as an MNP in my base VCF.

Is CountVariants running additional processes that I'm missing?

Thank you!

Issue · Github
by Geraldine_VdAuwera

Issue Number
1905
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
ronlevine

Best Answer

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Can you please post your CountVariants command lines, as well as the output you get?
  • JoniColemanJoniColeman LondonMember

    Command lines:

    java \
    -jar /humgen/gsa-hpprojects/GATK/bin/current/GenomeAnalysisTK.jar \
    -T VariantEval \
    --eval INPUT.vcf \
    -R /humgen/gsa-hpprojects/GATK/bundle/current/b37/human_g1k_v37.fasta \
    --dbsnp /humgen/gsa-hpprojects/GATK/bundle/current/b37/1000G_phase3_v4_20130502.sites.vcf \
    -nt 2 \
    -noEV \
    -noST \
    -ST Sample \
    -ST Filter \
    -ST Novelty \
    -EV CountVariants \
    -EV CompOverlap \
    -EV TiTvVariantEvaluator \
    -o OUTPUT_FILE
    

    Example output from CountVariants

    CountVariants   CompRod EvalRod Filter  Novelty Sample  nProcessedLoci  nCalledLoci nRefLoci    nVariantLoci    variantRate variantRatePerBp    nSNPs   nMNPs   nInsertions nDeletions  nComplex    nSymbolic   nMixed  nNoCalls    nHets   nHomRef nHomVar nSingletons nHomDerived heterozygosity  heterozygosityPerBp hetHomRatio indelRate   indelRatePerBp  insertionDeletionRatio
    CountVariants   dbsnp   eval    called  all 09C100188   3101804739  17716627    15823845    1892782 0.00061022  1638    1651467 3193    87258   97177   53074   0   613 33762   1153527 15790083    739255  10604   0   3.72E-04    2688    1.56    7.66E-05    13059   0.9
    CountVariants   dbsnp   eval    called  all 10C107613   3101804739  17716627    15834180    1882447 0.00060689  1647    1648755 3228    86741   94339   48770   0   614 35835   1142119 15798345    740328  8173    0   3.68E-04    2715    1.54    7.41E-05    13494   0.92
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Thanks, I think this might be a bug, but will consult the team to make sure. Would you be able to provide us with a snippet of your VCF that reproduces the issue?

  • JoniColemanJoniColeman LondonMember

    Thanks - attached is Example.zip. This contains:

    INPUT.recode.vcf, which is one of my samples on chr22;
    INPUT_Annotated.recode.vcf, which is that file annotated by VariantAnnotator;
    INPUT_Count_Var, which is the output from Count variants as run above.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @JoniColeman
    Hi,

    Thanks, I will have a look soon.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator
    edited May 2017

    @JoniColeman
    Hi,

    Sorry for the delay. I am not sure what exactly is going on here, but it looks like a bug in both tools. I will talk to the developers and see what they say. I may just be missing some extra flag. I should get back to you within the next week.

    -Sheila

    EDIT: It looks like VariantEval's CountVariants is doing the correct thing. It is counting MNPs only if the phased variants are all present in the interval. However, I will still need to check with the developers about why SelectVariants is not outputting the sites in its output.

  • JoniColemanJoniColeman LondonMember

    Thank you, Ron, that's very helpful - I hadn't appreciated the all alleles versus individual sample distinction in the tools.

    A couple of comments:

    AAAGAAAGAAAG <-> GAAGAAAGAAAG appears to only change at a single nucleotide (i.e. it could be written A <-> G if you weren't also listing the indel at the same locus) - is this considered an MNP rather than a SNP?

    SelectVariants doesn't seem very useful for this purpose - for example, in my data there is presumably no site where there is only an MNP (that is, SelectVariants must categorise all of them as MIXED of various forms, implying other variant types at all of these loci). That may be a true genomic feature - but it does make identifying these sites somewhat tricky!

    Thanks again!

  • ronlevineronlevine Lexington, MAMember
    edited June 2017

    @JoniColeman SelectVariants will trim the alleles to their common minimal representation. The spanning deletion (*) is preventing AAAGAAAGAAAG* GAAGAAAGAAAG from being represented as A* G.

Sign In or Register to comment.