On Monday and Tuesday, November 12-13, the communications team will be out of the office for a U.S. federal holiday and a team event. We will be back in action on November 14th and apologize for any inconvenience this may cause. Thank you for using the forum.

VQSR and VariantAnnotator on Samtools VCFs

Hi everyone!
My goal is to run VQSR on VCFs generated with samtools mpileup.
According to GATK best practices first i have to run VariantAnnotator on each of my VCFs in order to do that.

here's the annotation options I included in the command line:

--annotation QualByDepth
--annotation RMSMappingQuality
--annotation MappingQualityRankSumTest
--annotation ReadPosRankSumTest
--annotation FisherStrand
--annotation StrandOddsRatio
--annotation DepthPerSampleHC
--annotation InbreedingCoeff

unfortunately it returns these warnings:

WARN 17:26:07,297 StrandBiasTest - No StrandBiasBySample annotation or read data was found. Strand bias annotations will not be output.
WARN 17:26:07,297 InbreedingCoeff - Annotation will not be calculated. InbreedingCoeff requires at least 10 unrelated samples.
WARN 17:26:07,297 StrandBiasTest - No StrandBiasBySample annotation or read data was found. Strand bias annotations will not be output.
WARN 17:29:25,076 AnnotationUtils - DP annotation will not be calculated, must be called from HaplotypeCaller or MuTect2, not VariantAnnotator

I don't care about DP as for WXS it is not required, but the others are mandatory for VQSR (right?)

I'm using samtools to generate VCF for a couple of reasons: i need to call variants from single samples, i need a tool that calls "everything" without filters in order to apply custom downstream filtering. the idea is to try first the VQSR approach for the artifact filtering instead of hardfiltering (which is more tricky and complex).
I was thinking about using Haplotype Caller (also because it would be easier to use HC VCFs with VQSR or other GATK tools), but for what I have understood it is meant to find "germline SNPs", which doesn't fit my needs as I am looking for novel SNV (mutations) and not SNPs (am I right?).

Thank you very much.

Issue · Github
by shlee

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

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @Puputnik, our VQSR instructions are formulated for variant calls produced by GATK (specifically HaplotypeCaller) and are NOT appropriate for filtering calls made by other callers.

    To be clear, what we consider "germline" variants includes de novo mutations that are present in all the cells of the individual. This covers mutations that would have occurred either in the germline cells / zygotes in the parent, or during the initial formation of the embryo. This is in contrast to "somatic" variants that arise at any later stage such that they are not uniformly present in all cells (like cancer mutations). So HaplotypeCaller is appropriate for your work, unless I'm misunderstanding what you want.

  • PuputnikPuputnik ItalyMember

    Hi @Geraldine_VdAuwera , thanks for the response.
    Unfortunately I'm interested in any kind of variant (including "late" cancerous mutation) irrespective of allelic frequency etc.
    In other words i need everything that differ from reference genome, in order to apply custom downstream filtering later.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Then you should really consider using a combination of callers: HaplotypeCaller for germline and Mutect2 for somatic variation. Each provides a modeling algorithm that is appropriate for the relevant class of variation. They are both very sensitive and allow you to apply further custom filtering.

  • PuputnikPuputnik ItalyMember

    Hi Geraldine!
    I'm considering it, however i'm worried about the fact that Mutect2 can be used only with paired samples (tumor and matched normal), while I need to launch it on single samples.
    I was wondering: what if i launch mutect in "artifact detection mode" on a single sample? it should return any variant (even artifacts, but i don't care as i will filter them later) right?.
    Also, it is possible to use VQSR on VCFs generated with mutect in this way? i have a couple of these VCFs that I used in the past to generate a PON and i saw that the annotations required for VQSR are missing, is it possible to include them with the command "--annotation" in mutect command line?

  • shleeshlee CambridgeMember, Administrator, Broadie, Moderator admin
    edited November 2017

    Hi @Puputnik,

    The GATK4 Mutect2 documentation (for vbeta.6 here) gives an example command for running a tumor sample sans matched normal. It is essentially the same as the single normal sample calling command for panel of normals creation. This is not something we recommend in somatic analyses, rather it is a feature of the tool that you could exploit for your particular aims.

    GATK4 Mutect2 differs from GATK3 MuTect2. GATK4 no longer requires the artifact detection mode parameter. According to current documentation, it seems you can add additional annotations, but I believe these are limited to those that make sense for Mutect2. Check out Mutect2's --annotation, --annotationGroup and --annotationsToExclude options.

    However, as Geraldine says, VQSR is inappropriate for calls made by other callers than HaplotypeCaller (typos corrected). So any forays in this direction would be unsupported by us.

    You should also know that GATK4 separates out the filtering functionality into a tool called FilterMutectCalls. This allows iterative tweaking of per-sample filtering parameters, based on annotations available to Mutect2. Perhaps these may be sufficient for your per-sample needs.

    Perhaps our recent Mutect2 hands-on tutorial (Somatic SNVs and Indels) would be helpful to you. The extended version I wrote is in the 1709 tutorial folder within the workshop materials and gives an example use of --annotation. A link to the latter is at https://software.broadinstitute.org/gatk/documentation/presentations.

    Post edited by shlee on
  • PuputnikPuputnik ItalyMember

    Hi @shlee,

    thank you for the suggestions!
    I'm approaching with Mutect, however I have a couple of problems:
    1) I see the new GATK beta has been released, however the download link is still providing the 4.beta.5 version.
    2) I see Mutect2 doesn't require the COSMIC database anymore, and now the "af-only-gnomad.vcf.gz" is required.
    For what I've understood they are used for the same purpose: to provide known variants for labeling purposes (in other words, both of them are whitelists). but what are the differences in using COSMIC instead of the gnomad dataset?
    Also, I've seen that the "af-only-gnomad.vcf.gz" is based on GRCh37/hg19, and a GRCh38 version isn't available yet.
    I need to run my analysis using GRCh38, so I can't use the "af-only-gnomad.vcf.gz", is it possible to run mutect2 without it? how will it affect the analysis?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Puputnik
    Hi,

    1) Yes, that should be fixed soon. In the meantime, you can download beta 6 here.

    2) af-only-gnomad.vcf.gz is actually to filter out possible germline variants, not to whitelist possible somatic mutations. You can read more about it in the tutorial Soo Hee pointed to above.

    You can liftover the b37 VCF to hg38 using Picard's LiftoverVcf.

    -Sheila

  • PuputnikPuputnik ItalyMember

    hi @Sheila

    I don't love Liftover very much as I've read that in some cases it can introduce errors while converting to hg38, however if it is the only way I will do so.

    Anyway, in the meantime I tried to launch Mutect2 on a single tumor sample without the gnomad file, and I see that every variant is annotated as "germline_risk", so there aren't any "PASS" variants in my VCF, is it normal? any suggestion?

    Thank you very much.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Puputnik
    Hi,

    Yes, that is expected if you do not input a germline_resource file. Have a look at this thread.

    -Sheila

  • Hello, As i have been reading above, I have RNAseq data on Human Hepatocytes. So to get a VCF files for my samples (GSM******), I should not use HaplotyCaller and use Mutect2 isntead?
    Also, How can I add the gene information (gene symbol, gene name) in the vcf generated from above?

    Regards,
    Anurag

  • KatieKatie United StatesMember

    Hi,
    I am just now switching to GATK 4.0. Is VariantAnnotator no longer supported? Thank you

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited January 12

    @anuragpassi
    Hi Anurag,

    I really depends what your end goal is. Can you tell us more about what you hope to accomplish?

    To add the annotations, you can use VariantAnnotator. Have a look at this thread.

    -Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Katie
    Hi,

    There are plans to port VariantAnnotator to GATK4. It just has not happened yet.

    -Sheila

    Issue · Github
    by Sheila

    Issue Number
    2862
    State
    open
    Last Updated
  • KatieKatie United StatesMember

    Thank you! Would you be able to let us know when that change is made?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Katie
    Hi,

    Yes, I will make a note to let you know when that happens :smile:

    -Sheila

Sign In or Register to comment.