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

Convert CombineVariant output to individual vcf files

HI,

I was curious if it possible to convert the output of CombineVariants back into individual sample vcf files. Based on answered questions in other threads, I am under the impression that you should not try to use VariantFiltration on a combined output if variants were not called on an previously merged file. Is that correct? If so, could I basically reverse CombineVariants after using SelectVariants to analyze samples individually?

Also, I would like to point out that that all the links to the "Tools documentation index" are broken and displaying 404. I am not if you are aware.

Thank you

Best Answers

Answers

  • Thank you Sheila!

    I was hoping to find a subset of SNPs that exist in some samples but not in others and then apply hard filtering to narrow down my list.

    I guess my problem really is that I am trying to look for somatic mutations using RNAseq (whole transcriptome) among clonal replicates of a non-model plant. After following the best-practices for variant calling for RNAseq, I am getting a lot of false negatives. Basically, HaplotypeCaller is saying that the alternate allele was not called when it clearly exists in the bam file when visualized in IGV. Because expression levels of RNA can vary greatly, sensitivity is crucial. I know now that my question above was not the best approach to sifting through a "rough draft" of SNPs (but I was under an quickly approaching deadline to get something out).

    I seem to keep backing myself into a corner every time I try something because it appears that my ultimate goal is a bit unique and causing me a lot of confusion. I know that some of the tools have not been adequately tested for RNAseq but what do you think the best practice would be to look for mutations among clones?

  • I have clonal meristem samples from one treatment group and one control group of a diploid perennial plant. We are trying to understand mutation accumulation under high stress in genes/regions that are susceptible to selective pressures. This is not related to tumors or cancer. Theoretically, all 12 samples should be the same genetically because they have been replicated via vegetative growth but have been subjected to different selective pressures.

    I have a reference but it is a non-model plant (i.e. not Arabidopsis). I have used bootstrapping methods for things like base recalibration.

    It had appeared to me that MuTech2 was fairly specific to cancer/tumors or humans. Is this true? Or is there a way to apply this pipeline in a different way? If not, is there a better way to implement other tools for this purpose?

    Again thank you for all your help!

  • shleeshlee CambridgeMember, Broadie, Moderator

    Hi @JmeAlena,

    Mutect2 workflows we provide are fairly specific to human cancer. However, the tool itself is versatile and should allow you to call differences between samples. We've just made public last Friday the GATK4 Mutect2 tutorial:

    You will not want to follow this workflow. I point it out because it describes Mutect2 tool features in detail in sections one and two. One way you can use Mutect2 is in tumor-only mode and by providing the control sample's calls as the --germline-resource. I'm in the process of writing up a blogpost to describe such a workflow. Keep an eye out on our blogs if you are interested.

  • That sounds great. Thanks for your advice! What is the --germline-resource?

  • Ok great! Thank you both for being so thorough in answering all of my questions. I look forward to your blog post.

    How would this be applied to comparing 12 samples to each other? Should they be previously merged or combined after? Or is this just used to satisfy the file requirements?

  • JmeAlenaJmeAlena Member
    edited January 31

    Also should this vcf with allele frequencies be called with HC before inputting to MuTech2? Because I have had some issues getting an empirical allele frequency in a working vcf. I have had the same issue as this forum thread.

    https://gatkforums.broadinstitute.org/gatk/discussion/comment/36428

  • shleeshlee CambridgeMember, Broadie, Moderator

    @JmeAlena,

    I'm rereading this thread and it seems you want to compare two groups of multiple samples to each other, not just individual samples to each other as I previously assumed.

    I have clonal meristem samples from one treatment group and one control group of a diploid perennial plant.

    Also, you say that

    Theoretically, all 12 samples should be the same genetically because they have been replicated via vegetative growth but have been subjected to different selective pressures.

    So your set up would be:

    • Control group: sample1, ...sampleN
    • Test group: sampleO, sampleP, sampleQ...sampleZ

    Presumably, you are looking for differences in expression and given you are using GATK, you are interested in differences in expressed alleles between your two groups. Is this the case?

    If you follow the Mutect2 tutorial sans a matched normal BAM, then for the --germline-resource you can provide a callset representing the control group. This germline resource can use actual allele frequencies representing the control group samples. You would analyze each of the test group's samples against the germline resource in independent Mutect2 runs.

    The caveat of this approach is that you will only call variant sites that are new to a sample (not present in the germline resource). Also, this approach doesn't call reversions to reference, which may be a change in the sample from the control group but which isn't considered somatic in the traditional sense.

    The yet-to-be-posted blog post outlines a way to compare one sample against another sample and call the contrasting alleles including reversions to reference. The blog post workflow is an off-label workflow--unofficial and unvetted. You should consider whether you want to somehow pool samples for each group or decide on a representative sample towards this alternate approach. If you decide to pool, then one approach is to do this after the first-pass tumor-only mode calling of each sample, using CombineVariants.

    I'm a bit concerned you say:

    I am getting a lot of false negatives. Basically, HaplotypeCaller is saying that the alternate allele was not called when it clearly exists in the bam file when visualized in IGV.

    Please make sure you can call the variants you expect with Mutect2. Both HaplotypeCaller and Mutect2 share code for haplotype assembly. As Sheila mentioned, if it the case that HaplotypeCaller isn't calling a low fraction allele, then switching to Mutect2 should help. One thing to be careful of when using assembly callers is intervals lists that have small intervals. If your setup doesn't allow you to call the expected variant, and you are applying an intervals list, I recommend trying calling with larger intervals, e.g. unsplit contigs, or no intervals. You can always subset out the variants in regions of interest post-calling.

  • Hi Shlee,

    I do want to compare samples to samples because the clones that I am using are all capable of mutating separately from each other and at low allele frequencies. I would like to be able to see called variants side by side in all samples to analyze in what samples they do and don't occur (similar to the viewing of combine variants). I do have samples that were not exposed to a harsh environment but mutation accumulation should not be ruled out in these samples. I am sorry, if my statements were misleading. Your original explanations seem like they may be effective to accomplish what I want.

    With all that said, can I use combine variants on my outputs from MuTech2?

    And I was curious that when I provide MuTech2 a vcf to call variants in my samples, how should I generate that vcf to provide empirical AF?

    Again, thank you for your thorough help and quick responses!

  • shleeshlee CambridgeMember, Broadie, Moderator

    Hi @JmeAlena,

    Yes, I have used CombineVariants to combine two Mutect2 callsets. However, be sure to perform all filtering, i.e. with FilterMutectCalls, prior to the combining. The combination is artificial so to speak so be wary of any annotations the tool may generate during the combination. Also, if two variants at the same site have different filtering status, then the default settings for CombineVariants will turn the site into a PASS site, which you may find undesirable, pending your use case for the combined VCF.

    Best to test what the tool does and do your own sanity checks as to whether the results meet expectations.

    The blogpost outlines steps for your 2nd question: https://software.broadinstitute.org/gatk/blog?id=11315.

  • Yes of course! Thank you! I am going to check this link out

Sign In or Register to comment.