How to force MuTect2 genotype all sites within intervals?

I cannot generate genotypes for all input sites with MuTect2 (-L input.vcf). I tried --output_mode EMIT_ALL_SITES and -gt_mode GENOTYPE_GIVEN_ALLELES without luck. Is there an equivalent option to --force_output in MuTect? I know it is more complicated with indels where it is not applicable to generate all possible indels, but I am only interested in the specific genotypes I used with -L.

Tagged:

Answers

  • kmhernankmhernan Chicago, ILMember

    I know some people have done this, but as far as I can tell they have their own GATK build. I have tried many many combinations of options in MuTect2 to do this without success. I definitely would like to do this especially to get the OxoG counts in variants from other somatic callers.

  • amjaddamjadd FinlandMember

    so if there is no way to do this, is there any plans to allow this functionality in future versions?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    No plans to do this at this time, but if more people tell us they want this functionality, we'll take that under consideration.
  • amjaddamjadd FinlandMember

    OK thanks Geraldine. I think it is an important feature. My use case is that I have multiple samples from the same individual so having reference and alternative read counts in each of the samples for each variant is crucial.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Ah, I see. I think there might be some work planned to be able to process multiple tumor samples from the same individual together, which might be a better way to deal with this question. I'll try to find out what is on the somatic development roadmap.
  • kmhernankmhernan Chicago, ILMember

    My use case is pretty much the same at the TCGA MC3... to get the OxoG data and some other statistics across multiple callers (e.g., varscan2, MuSE, somaticsniper, etc.) so that more standardized filtering can be performed downstream.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Some of that might be annotatable through VariantAnnotator -- we'll certainly do what we can to make that work. For the rest, I hope you'll understand it's difficult for us to prioritize feature requests that are motivated by purely extrinsic needs. Not that we don't see the value of such approaches -- we definitely support the goal of standardizing methods -- but there is already so much to do within the MuTect-centric pipeline, it's a hard sell to allocate effort to these peripheral goals.
  • alexcharneyalexcharney Mount SinaiMember

    Also tried many combinations of options to MuTect2 to generate OxoG metrics for variants called by other callers (e.g., MuTect1), but to no avail. If anyone has any luck or update on this front would be much appreciated.

  • alongaloralongalor Member
    edited April 9

    Also have a similar use case - have multiple samples for the same individual and would like to have "reference and alternative read counts in each of the samples for each variant" as stated earlier, in addition to other metrics. Is this possible with -gt_mode in Mutect2 GATK4? Or are there any other flags you would suggest using to glean any further insight about why a variant is called in one sample but not another?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alongalor
    Hi,

    There is now an option to emit germline sites in Mutect2 in the latest release. https://github.com/broadinstitute/gatk/pull/4522

    There is also an option to have GGA mode work in Mutect2 that is in review and should be in the next release. https://github.com/broadinstitute/gatk/issues/4555

    -Sheila

  • alongaloralongalor Member
    edited April 10

    Hi Sheila,

    Thanks a lot for these great tips! Will be sure to check both of these out (the latter as it becomes available!)

    I also noticed the --forceActive option and used it in conjunction with the -L argument, and found the output very useful!

    A question regarding to that: Say I am interested in a variant that occurs at position 1:15094 that Mutect2 did not emit for a particular sample.

    I found that by using the --forceActive and -L 1:14500-15500, the resulting VCF file emmits the variant at position 1:15094 (which is fantastic!).

    I was wondering what would be the narrowest band you recommend setting the --L arugment to in order to capture the variant, while maintaining the integrity of the "FILTER" column (I am hoping to reduce runtime for a large list of variants, while maintaining the validity of the "FILTER" reason?

    It would seem that give the following information:

    Maximum region size: 300 bp
    Extension increments: 100 bp
    

    -L chrom:pos_of_interest-400:pos_of_interest+400 might be prudent? Unless it is possible to have multiple extension increments? Do you know if this is possible?

    Thanks a lot,

    Alon

    Post edited by alongalor on
  • Update: Turns out this doesn't work well (if at all). I would also be grateful if the feature discussed above (some sort of --force_output) option would be available.

  • Update: Looks like the originally described request works in Mutect2 GATK4! Thanks a lot!

  • shleeshlee CambridgeMember, Broadie, Moderator

    Thanks for the clarification that it works well in GATK4 Mutect2 @alongalor.

  • Hi @shlee,

    After many tries to try to get this working, I was unfortunately not able to get all the sites in my vcf file passed in the --alleles flag to show up in the Mutect2 output. There is some overlap between the two, however I would like the overlap to be 100%. As you can see in my full command below, I also use the --genotyping-mode GENOTYPE_GIVEN_ALLELES and --output-mode EMIT_ALL_SITES flags.

    /n/data1/hms/dbmi/park/alon/software/gatk/gatk-4.0.3.0/gatk Mutect2 \
         -R /n/data1/hms/dbmi/park/alon/software/gatk/human_g1k_v37_decoy/human_g1k_v37_decoy.fasta \
         -I /n/data1/hms/dbmi/park/DATA/BSMCommonExperiment/ReferenceTissueProject/sourceData/sourceDataNew/bam/Walsh-Park/LBRCE-pfc-1b123.bam \
         -tumor LBRCE-pfc-1b123 \
         -I /n/data1/hms/dbmi/park/DATA/BSMCommonExperiment/ReferenceTissueProject/sourceData/sourceDataNew/bam/Chess-Walsh-Akbarian/Fibro_Common_7_DO16090243-final-all.bam \
         -normal Fibro_Common_7 \
         --germline-resource /n/data1/hms/dbmi/park/alon/software/gatk/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz \
         --genotyping-mode GENOTYPE_GIVEN_ALLELES \
         --alleles /n/data1/hms/dbmi/park/alon/files/BSMNCE/Mutect2Tests/Sestan_PASS.recode.150.vcf \
         --output-mode EMIT_ALL_SITES \
         --native-pair-hmm-threads 20 \
         -L /n/data1/hms/dbmi/park/alon/files/BSMNCE/Mutect2Tests/Sestan_PASS.recode.150.vcf \
         -O /n/data1/hms/dbmi/park/alon/files/BSMNCE/Mutect2Tests/LBRCE-pfc-1b123_v_Fibro_Common_7_DO16090243-final-all_GGA.vcf
    
    

    The vcf I pass to the --alleles argument is a Mutect2 vcf file of a bulk WGS sequenced tissue (a biological replicate of the sample resulting in the input bam I pass into the below command)

    As stated earlier, my use case is as follows: I have multiple samples for the same individual and would like to have "reference and alternative read counts in each of the samples for each variant", in addition to other metrics, in order to understand why some variants were called in one sample but not another.

    I would appreciate any help in solving this issue. As far as I understand from the documentation, the -O file should contain sites at all the positions passed into the -alleles argument.

    Thanks,

    Alon

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alongalor
    Hi Alon,

    I think in your original post, you are asking about -ip.

    -Sheila

  • alongaloralongalor Member
    edited April 14

    Hi Sheila, Thanks for your response. That is certainly helpful, however the main issue at hand here is that I cannot get Mutect2 to output a line in the output vcf file for every variant i pass in the -alleles argument. I have found that adding the below flags to my above posted command increase the percentage of lines in the output vcf file emitted for every variant i pass in the -alleles argument from ~20% to ~80%, however I really do need 100% coverage and am wondering if there are any additional flags I can use to make that happen!

    --standard-min-confidence-threshold-for-calling 0 \
    --tumor-lod-to-emit -100000000000 \
    --min-base-quality-score 0 \
    --initial-tumor-lod -100000000000 \
    

    Also I should add that optimally I think it is probably best not to use these flags if possible since they could potentially impact the "FILTER" reason.

    Thanks a lot,

    Alon

  • SheilaSheila Broad InstituteMember, Broadie, Moderator
  • Hi Sheila, thanks a lot for the input - unfortunately I still do not get ~20% of sites. This is really important for our analysis so it would be amazing if it would be possible get get 100% of the sites passed in the -alleles arg! thanks a lot for your help!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alongalor
    Hi Alon,

    Okay. I have two more arguments for you to add :smile:

    --genotype-germline-sites and --genotype-pon-sites will help with sites that may be in the PoN or possible germline events.

    I hope those do the trick!

    -Sheila

  • Hi Sheila,

    Thanks so much for your continued help with this! These are all great ideas! Unfortunately these still don’t get me past 85% :( would you have any additional filters to recommend or do you think there might be a flag we are somehow missing that will emit all sites passed to it like a combination of --genotyping-mode GENOTYPE_GIVEN_ALLELES and --output-mode EMIT_ALL_SITES -alleles regions.vcf is supposed to provide, as explained by the documentation.

    Thanks again for your help with this,

    Alon

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alongalor
    Hi Alon,

    To confirm, the issue is that when you pass in a VCF in GGA mode, you are missing ~15% of sites in the final VCF? If so, can you post some sites in the GGA VCF that are missing in the output VCF?

    Thanks,
    Sheila

  • alongaloralongalor Member

    Hi Sheila,

    That is correct.

    When I provide a bed file of 10 sites:

    1       15093   15094
    1       82114   82115
    1       102333  102334
    1       241193  241194
    1       568360  568361
    1       646800  646801
    1       702943  702944
    1       1595952 1595953
    1       2581331 2581332
    1       2585648 2585649
    

    The following 2 are missed in the output vcf:

    1       646800  646801
    1       1595952 1595953
    

    Thanks again for your help!

    Alon

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alongalor
    Hi Alon,

    Can you post the BAM file and bamout file for those two missed sites? Can you also post the BAM and bamout for 1 site that is called?

    Thanks,
    Sheila

  • alongaloralongalor Member

    Hi Sheila,

    Unfortunately the BAM file is huge and so I wont be able to. However, I can post the bamout file if that would be sufficient.

    Thanks for your help,

    Alon

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @alongalor
    Hi Alon,

    I didn't mean the entire BAM files :smile: I just meant IGV screenshots. I would like to see ~300 bases before and after the sites of interest to see if there is anything funny going on there.

    Thanks,
    Sheila

  • cgibsoncgibson Member

    Just wondering if there had been any update on this thread...I also need to pass a list of alleles to M2 for evaluation (for the same reason as some others on this thread, need to get read pair orientation data for all variants called by a different variant caller) and was about to try this strategy...but if it doesn't work for all variants I may try something else. Thanks!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @cgibson
    Hi,

    No, I don't think there was any progression after my last post. Any chance you can try out all the recommendations and report back what you find?

    Thanks,
    Sheila

Sign In or Register to comment.