pon_filter get more sites? different bed file?

hi, I get more sites after pon_filter.
in the line of the more sites in vcf, there is a label "panel_of_normals", and the big bed contains all the regions of small bed,
I tried find a detailed interpret of somatic vcf, but just find a old version, is there a new version for gatk4?
thanks a lot.

Tagged:

Best Answer

  • manbamanba ✭✭
    edited December 2018 Accepted Answer

    Q16: My most puzzle is why after pon filter, there are some new variant site appearing?
    the following is an example
    (chr7 128845018 . C T . germline_risk DP=148;ECNT=1;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;TLOD=260.68 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:OBAM:OBAMRC:OBF:OBP:OBQ:OBQRC:SA_MAP_AF:SA_POST_PROB 0/1:69,79:0.534:69,79:0,0:34:177,177:60:20:false:false:.:.:36.73:100.00:0.00,0.535,0.534:0.029,0.025,0.946)

    and after pon filter(I mean the vcf still labeled with germline_risk. panel of normals ....), there are many files really decrease the variant site number, I am curious about why these dropped, but most sites labled, really amazing. the follwoing is an example of the sites dropped after pon
    (
    chr4 55602765 . G A,C . germline_risk DP=314;ECNT=1;POP_AF=1.000e-03,1.000e-03;P_GERMLINE=-2.807e-03,-2.746e-04;TLOD=4.92,6.57 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:OBAM:OBAMRC:OBF:OBP:OBQ:OBQRC:SA_MAP_AF:SA_POST_PROB 0/1/2:301,4,5:0.020,0.023:300,4,5:0,0,0:34,36:197,197,197:60,60:38,46:false:false:.,.:.,.:100.00:36.73:0.020,0.00,0.016:3.071e-03,5.067e-03,0.992
    )

    to be more clear, I give a screenshot of before and after Pon.

    so I want to know after PoN. which kind of sites will be labeled with words like "panel_of_pon", which kind of site will be dropped, which kind of site will new generation.

    in this https://software.broadinstitute.org/gatk/blog?id=10911, it said

    "
    GATK3 MuTect2 prefilters sites in the germline resource regardless of the allele in the tumor. GATK4 Mutect2 distinguishes alleles in the germline resource and only filters the site if the tumor allele matches. If the alleles are different, then the tool considers the allele a putative somatic mutation.

    Filtering of sites in the panel of normals (PoN) and the matched normal remains unchanged, except that the tool will prefilter most of these such that site records are absent from the VCF.
    "
    so I think lines containing words like "panel_of_pon" should be filtered manually, but if just as your standard, manually dropped if contaning any words like "clustered_events;germline_risk;read_position;t_lod", there will be no sites, really upset.

Answers

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    @manba GATK 4 Mutect2 emits regular spec-compliant VCFs. I'm not sure what a somatic VCF is. Also, in GATK 4 the panel of normals is a VCF, not a bed.

    I'm afraid none of us can be very helpful for questions involving previous versions of Mutect2. We have validated the superiority of GATK 4 Mutect2 very thoroughly and since it's such an obvious difference we haven't bothered to familiarize ourselves with old versions.

  • manbamanba Member ✭✭
    edited December 2018

    I used the gatk4.0.0.0.

    a somatic vcf means the vcf is called by mutect2 in gatk4 after some filters(sorry for the not clearly definition ). Different bed files means I used a bed file contains more genes to create pon for
    a lot of bed files(these contains less genes )

    let me raise my puzzle by ordered questions, thanks a lot.
    I mainly focus on website https://gatkforums.broadinstitute.org/gatk/discussion/11136/how-to-call-somatic-mutations-using-gatk4-mutect2 and GATKwr23-S1-Somatic_SNVs_and_indels.pdf in your workshop.

    Q1: I am puzzled about one sentense in the following picture, "All 17 of these happen to be filtered.", how do you understand filtered, means not appear in the vcf after pon filter or just add some label. and in your workshop pdf, you listed some filters, does it means if the line appear any of the 15 rules(FilterMutectCalls filters for multiple criteria), these vcf sites should be dropped manually, I do not know the filtered in your vcf is manual or automatical , also, it happened to see all the 17 filtered sites in your vcf meet some of the 15 15 rules(FilterMutectCalls filters for multiple criteria).

    Q2: Through https://software.broadinstitute.org/gatk/blog?id=10911, I know "GATK4 Mutect2 distinguishes alleles in the germline resource and only filters the site if the tumor allele matches. If the alleles are different, then the tool considers the allele a putative somatic mutation. means it not simply filtered sites in pon, am I right, but I am surprised to find that some tumor vcf files filtered one site in pon, while some other not filter this site, and all the tumot vcf files in this site is the same (G -> A), so I really could not understand this.

    Q3: I also surprised to find that some vcf files after pon filter generate some new sites in vcf, really amazing.

    Q4: do you think is the error in gatk4.0.0.0, but can be fixed in the latest version?

    Q5: can you give some advice in filtering variant in mutect2, for example, to lower or higher some parameters in https://software.broadinstitute.org/gatk/documentation/tooldocs/4.0.6.0/org_broadinstitute_hellbender_tools_walkers_mutect_FilterMutectCalls.php

    --heterozygosity
    --heterozygosity-stdev
    --indel-heterozygosity
    --interval-merging-rule
    --log-somatic-prior
    --max-alt-allele-count
    --max-contamination-probability
    --max-events-in-region
    --max-germline-posterior
    --max-median-fragment-length-difference
    --max-strand-artifact-probability
    --min-base-quality-score
    --min-median-base-quality
    --min-median-mapping-quality
    --min-median-read-position
    --min-pcr-slippage-size
    --min-strand-artifact-allele-fraction
    --native-pair-hmm-threads
    --native-pair-hmm-use-double-precision
    --normal-artifact-lod
    --normal-p-value-threshold
    --num-reference-samples-if-no-call

    Q6: Why some in 15 rules has no value? like germline_risk, t_lod, no exact value , like 0.1, 0.2, thanks a lot

  • manbamanba Member ✭✭

    Q6: I used gatk4.0.0.0 all the time, but the $7 in my twicefiltered.vcf is not word like "PASS', it like the following

    germline_risk
    base_quality;germline_risk;read_position;strand_artifact;t_lod
    germline_risk;multiallelic
    germline_risk
    germline_risk
    germline_risk

    and in the following picture, Of the filtered calls, ~56% (1,694/3,022) are filtered singly. if I not incorrectly understand, the third column solo filter instances means sites that are dropped, so if o appear in that column, means it filter no sites? am I right?

    about the author's question "Which filters appear to have the greatest impact? What types of calls do you think compels manual review?
    "
    from the aspect of absolute value, 'panel of normals' has the greatest impact,
    from the aspect of ratio, orientention_bias (100% dropped) has the greatest impact.
    am I right? so why the sample labels, some dropped, but some lucky kept, is there some backgroud reasons?

    the author' s question " What types of calls do you think compels manual review?" is also what I want to ask. really a intesting question, can you tell me the answer of your advice.

    can not be more thanksed, Thank you so much.

  • manbamanba Member ✭✭

    Q7: how can I know each filter filter how many sites just like the above picture? thanks a lot

  • manbamanba Member ✭✭
    edited December 2018

    https://github.com/broadinstitute/gatk/blob/master/docs/mutect/mutect.pdf


    8: A more surprised discovery: In the above link about concrete methods of Mutect2(it signed data 2018.11.20)
    it said remove all panel of normals label, no matter in https://gatkforums.broadinstitute.org/gatk/discussion/11136/how-to-call-somatic-mutations-using-gatk4-mutect2 or my results, I never find remove all.

    urgent for help. thanks a lot

  • manbamanba Member ✭✭

    Q10: the most important question is after FilterMutectCalls, the label in each line is a result of after dropping some varaint site or not (I think is yes), not solving this problem, other questions can not be solved too

  • manbamanba Member ✭✭

    Q11: Filtering thresholds for both normal-artifact-lod (default threshold 0.0) and tumor-lod (default threshold 5.3) can be set in this tool(FilterMutectCalls). the value is a log-10 value, so it means 1 reads support normal, 10^5.3 reads support homref hypothesis(according to website https://gatkforums.broadinstitute.org/gatk/discussion/comment/54688#Comment_54688)

  • manbamanba Member ✭✭

    Q12: Still the parameter of FilterMutectCalls, should we set the mode? if select the EMIT_ALL_CONFIDENT_SITES, how is the confident defined? and it said "produces calls at variant sites and confident reference sites
    "
    means only confident reference sites, but not confident variant sites ? thanks a lot

  • manbamanba Member ✭✭

    Q12: is it possible to ignore some sites to be filtered in Filtermutectcalls? because it is hotpot, I want to keep. thanks a lot

  • manbamanba Member ✭✭

    Q13: form DP=152;ECNT=1;IN_PON;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;TLOD=484.72,
    the TLOD is bigger than the default 5.3, so it pass the filter,

    and IN_PON means this site in pon vcf,

    ECNT means cluster_event(max-events-in-region is the maximum allowable number of called variants co-occurring in a single assembly region. If the number of called variants exceeds this they will all be filtered.) so here is 1, not exceeds, but what is cutoff

    DP=152;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04 cutoff is what, especilly what they stand for in your mehods about pdf.

  • manbamanba Member ✭✭

    Q14: The filter column can not find PASS, means there is no site can be treated as true variant, I looked the vcf before pon and after pon(after Filtermutect2call command), both has no PASS in Format column, does it really means there is no true vaiant site,

    urgent , thanks a lot

  • manbamanba Member ✭✭
    edited December 2018

    Q15 I am sorry for misunderstanding this picture. both column 2 and 3 stands for filter sites number.

    so as your standard, if any filter rules appear in the line, the site is not real? or we should combine the DP and AF to compare with pon(if AF is much difference, considered to be a somatic varainr)

    Do you think is too strict to just select the value of column 'FILTER' is 'PASS' as true variant?
    '

  • manbamanba Member ✭✭
    edited December 2018 Accepted Answer

    Q16: My most puzzle is why after pon filter, there are some new variant site appearing?
    the following is an example
    (chr7 128845018 . C T . germline_risk DP=148;ECNT=1;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;TLOD=260.68 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:OBAM:OBAMRC:OBF:OBP:OBQ:OBQRC:SA_MAP_AF:SA_POST_PROB 0/1:69,79:0.534:69,79:0,0:34:177,177:60:20:false:false:.:.:36.73:100.00:0.00,0.535,0.534:0.029,0.025,0.946)

    and after pon filter(I mean the vcf still labeled with germline_risk. panel of normals ....), there are many files really decrease the variant site number, I am curious about why these dropped, but most sites labled, really amazing. the follwoing is an example of the sites dropped after pon
    (
    chr4 55602765 . G A,C . germline_risk DP=314;ECNT=1;POP_AF=1.000e-03,1.000e-03;P_GERMLINE=-2.807e-03,-2.746e-04;TLOD=4.92,6.57 GT:AD:AF:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:OBAM:OBAMRC:OBF:OBP:OBQ:OBQRC:SA_MAP_AF:SA_POST_PROB 0/1/2:301,4,5:0.020,0.023:300,4,5:0,0,0:34,36:197,197,197:60,60:38,46:false:false:.,.:.,.:100.00:36.73:0.020,0.00,0.016:3.071e-03,5.067e-03,0.992
    )

    to be more clear, I give a screenshot of before and after Pon.

    so I want to know after PoN. which kind of sites will be labeled with words like "panel_of_pon", which kind of site will be dropped, which kind of site will new generation.

    in this https://software.broadinstitute.org/gatk/blog?id=10911, it said

    "
    GATK3 MuTect2 prefilters sites in the germline resource regardless of the allele in the tumor. GATK4 Mutect2 distinguishes alleles in the germline resource and only filters the site if the tumor allele matches. If the alleles are different, then the tool considers the allele a putative somatic mutation.

    Filtering of sites in the panel of normals (PoN) and the matched normal remains unchanged, except that the tool will prefilter most of these such that site records are absent from the VCF.
    "
    so I think lines containing words like "panel_of_pon" should be filtered manually, but if just as your standard, manually dropped if contaning any words like "clustered_events;germline_risk;read_position;t_lod", there will be no sites, really upset.

  • manbamanba Member ✭✭

    new sites in vcf after pon filter means not present in the PoN.vcf made by many normals(the eight column vcf file, I think you know my meaning) and not present in the vcf after Filtermutect2calls before pon filter.

    really strange.

  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018

    now I find most default values of the following criterion of the above picture, but find some unclear events.

    the following 6 criterion seems to have no value in the column of INFO even the key appear in the column of FILTER

    _strand_artifact
    base_quality
    mapping_quality
    read_position
    panel_of_normals
    str_contraction _

    base_quality;clustered_events;germline_risk;read_position;t_lod DP=83;ECNT=3;POP_AF=1.000e-03;P_GERMLINE=-5.804e-02;TLOD=3.54

    germline_risk;str_contraction DP=3394;ECNT=2;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;RPA=4,3;RU=CTGT;STR;TLOD=53.01

    germline_risk;panel_of_normals DP=128;ECNT=1;IN_PON;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;TLOD=22.27

    clustered_events;germline_risk;read_position DP=73;ECNT=3;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;TLOD=10.15

    clustered_events;germline_risk;mapping_quality;panel_of_normals DP=96;ECNT=5;IN_PON;POP_AF=1.000e-03;P_GERMLINE=-2.169e-04;TLOD=13.74

    base_quality;germline_risk;read_position;strand_artifact DP=131;ECNT=1;POP_AF=1.000e-03;P_GERMLINE=-2.404e-04;TLOD=6.96

    the criterion artifact in normal seems neither appear in the column of FILTER **or **INFO.

    the value of germline_risk is minus, but why a p value can be minus? and what _POP_AF _stands for , what the_ AF_ stands for in the column of FORMAT

    the most serious question is in https://software.broadinstitute.org/gatk/documentation/tooldocs/4.0.11.0/org_broadinstitute_hellbender_tools_walkers_mutect_FilterMutectCalls.php,
    you said exceeds or lower than the value of criterion, it will be filtered, but never mentioned just add a lable in the columnof FILTER, just like "base_quality;germline_risk;read_position;strand_artifact"

    see the red line , just said filter, in fact, not the thing

    gatk4 said it give the filter resposibilty to FilterMutectCalls, just not like a mix in gatk3.
    I want to ask is there a complete seperate from Mutect2 of the filter function?

    thanks a lot.


  • means all in pon but high fraction in matched normal, why, pon also made of normal samples, if I also use the matched normal in the pon, how will you do, your Mutect2 also support add pon in the paired mode? is it a bit inconsistent?

    thanks a lot

  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018


    if like you said, pon can not replace matched normal, if I put my matched normal into pon, how how how you treat.

    pon just for filter sequence error as the bold words said?
    we focus on human, why lack we a common germline variant resource? it supplied in genomeAD , and you also used that, why we lack, obviously we do not lack. and I think it really need a clear classfication compare between pon and a common germline variant resource. to give a better definition about what they really do each other. do not you think it has a conflict with your pdf in your workshop.

    means if there is not matched sample, the pon has no ability to filter gemline varaint ? and if really so, I think it is realy a urgent time for you to call up that genomeAD must be used.

    thanks a lot. due to my low ability, I must have misunderstood a lot of things, hopefully, all the staff active in gatk forum can forgive my ignorance and folly

  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018

    I may found another error in the doc. pictures comes from website
    https://gatkforums.broadinstitute.org/gatk/discussion/12986/af-of-alleles-not-in-resource-in-tumor-only-model-in-gatk4-mutect2

    https://software.broadinstitute.org/gatk/documentation/tooldocs/4.0.11.0/org_broadinstitute_hellbender_tools_walkers_mutect_Mutect2.php#--normal-lod. _ is this really default is a minus?_

    @shlee , since you mentioned this may change according to version. how can I get the real default value?

    and I think we just get the "--af-of-alleles-not-in-resource " from genomeAD website. or your broad googleclound. is the resoure in your broad googleclound made by yourself or download from genomeAD?

    is there a easy method for me to know how many samples involved in my "--af-of-alleles-not-in-resource " ?

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    @manba

    I used the gatk4.0.0.0.

    This version is one year out of date. Using the most recent release is almost always a good idea.

    Q1: I am puzzled about one sentense in the following picture, "All 17 of these happen to be filtered.", how do you understand filtered, means not appear in the vcf.

    Mutect2 output, like that of all GATK 4 tools, conforms to the VCF spec: https://samtools.github.io/hts-specs/VCFv4.2.pdf

    Q3: I also surprised to find that some vcf files after pon filter generate some new sites in vcf, really amazing.

    By default the we skip assembly and variant calling for sites in the pon, so a pon site might be absent from the called vcf.

    Q4: do you think is the error in gatk4.0.0.0, but can be fixed in the latest version?

    I do not think this behavior is an error.

    Q5: can you give some advice in filtering variant in mutect2

    Use the defaults unless you really, really, know what you are doing.

    Q6: Why some in 15 rules has no value? like germline_risk, t_lod

    I'm not sure what you mean because those filters have adjustable command line thresholds.

    Q7: how can I know each filter filter how many sites just like the above picture? thanks a lot

    You can find the number of sites uniquely filtered by, for example, clustered events, with grep -e "[^;]clustered_events[^;]" calls.vcf | wc -l

    8: A more surprised discovery: In the above link about concrete methods of Mutect2(it signed data 2018.11.20)

    it said remove all panel of normals label, no matter in https://gatkforums.broadinstitute.org/gatk/discussion/11136/how-to-call-somatic-mutations-using-gatk4-mutect2 or my results, I never find remove all.

    See answer to Q1.

    Q10: the most important question is after FilterMutectCalls, the label in each line is a result of after dropping some varaint site or not (I think is yes), not solving this problem, other questions can not be solved too

    See answer to Q1. FilterMutectCalls applies filters in the filter column in accordance with the VCF spec.

    Q11: Filtering thresholds for both normal-artifact-lod (default threshold 0.0) and tumor-lod (default threshold 5.3) can be set in this tool(FilterMutectCalls). the value is a log-10 value, so it means 1 reads support normal, 10^5.3 reads support homref hypothesis(according to website

    A likelihood is not the same as a count.

    Q12: Still the parameter of FilterMutectCalls, should we set the mode?

    See answer to Q5.

  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018

    @davidben realy thank you, sorry to disturb you in your holiday. but I really do not agree with what you said, I hope you can see all my answers in this post before you answers me, thanks a lot.

    It is always not easy to update to the latest because we need to keep the stability of the result at that period. we also need to be responsible for the patient.
    also I do not think the result is just the reason of old version.

    Q3: I also surprised to find that some vcf files after pon filter generate some new sites in vcf, really amazing.

    By default the we skip assembly and variant calling for sites in the pon, so a pon site might be absent from the called vcf.

    it is obvious what you said is not what I want to know. according to the doc by shlee, skip assembly and variant calling for sites in the pon, but in fact, what i said is after pon, there some new variant sites, new means " not" in the vcf before pon, but appear after pon.

    _and I think Filtermutect2calls is completele not automatically filter, because there are a lot of varaint labled with "t_lod, germline_risk,panel_of_normals", if really like what you said "skip assembly and variant calling for sites in the pon", I think it will not have the label, at least not "panel_of_normals".__**

    and if the matched normal is also in pon, how will you do?
    "Another major difference is in site versus allele filtering against the germline resource. GATK3 MuTect2 prefilters sites in the germline resource regardless of the allele in the tumor. GATK4 Mutect2 distinguishes alleles in the germline resource and only filters the site if the tumor allele matches. If the alleles are different, then the tool considers the allele a putative somatic mutation."

    we alos need to consider if the variant in germline is both different from reference.

    Q6: Why some in 15 rules has no value? like germline_risk, t_lod

    I'm not sure what you mean because those filters have adjustable command line thresholds.

    here I mean the column of FILTER and INFO, you can see, there are some label have no corresponding value in the INFO column.
    for example, can you see the value of "base_quality" in the INFO column, and do you know, what POP_AF correspond to? if you know, do you know why P_GERMLINE is minus, and what the difference between POP_AF and AF?

    base_quality;clustered_events;germline_risk;read_position;t_lod DP=83;ECNT=3;POP_AF=1.000e-03;P_GERMLINE=-5.804e-02;TLOD=3.54
    "

    if it is ok. hope we can have a talk in skype, thanks a lot. I hope there is no Misunderstanding between us. thanks a lot

    Post edited by picard_gatk_mj on
  • Q11: Filtering thresholds for both normal-artifact-lod (default threshold 0.0) and tumor-lod (default threshold 5.3) can be set in this tool(FilterMutectCalls). the value is a log-10 value, so it means 1 reads support normal, 10^5.3 reads support homref hypothesis(according to website

    A likelihood is not the same as a count.

    so can you help me finding out the relationship between likelihood and count in this. thanks a lot

  • manbamanba Member ✭✭

    The germline resource gatk offered, the_ chromosome name does not start with "chr"_, we need to convert it by myself. the following command is ok

    sed -ri '/^\s*#/!s#^#chr#' af-only-gnomad.raw.sites.vcf

    after I apply the germline resource, there is no variant anymore, means no matter the vcf with or without pon, just end with the column "#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample_name ".

    what does it imply, still no good variant, so I surprised to find the argument --germline-resource is too powerful, does it contain all the rules in the command of FilterMutectCalls . rules means the following picture.

    thanks a lot

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    It is always not easy to update to the latest because we need to keep the stability of the result at that period. we also need to be responsible for the patient.

    @manba The Mutect2 official workflow costs about $1 for a tumor-normal WGS pair. I would strongly recommend re-running old results with the most recent release. If you want, you could wait until the 4.1 release in January 2019, which will have several new advances, including a multi-sample mode.

    also I do not think the result is just the reason of old version.

    There have been a lot of improvements in the last year. Mutect2 is a mature tool, and we are very careful to evaluate it before every release, but it is also under very active development. 4.0.0.0 is now obsolete.

    it is obvious what you said is not what I want to know. according to the doc by shlee, skip assembly and variant calling for sites in the pon, but in fact, what i said is after pon, there some new variant sites, new means " not" in the vcf before pon, but appear after pon.

    Please supply an example along with your command lines used for Mutect2 and FilterMutectCalls.

    _and I think Filtermutect2calls is completele not automatically filter, because there are a lot of varaint labled with "t_lod, germline_risk,panel_of_normals", if really like what you said "skip assembly and variant calling for sites in the pon", I think it will not have the label, at least not "panel_of_normals".__**

    Sometimes an assembly region triggered by a site that is not in the panel of normals will include a panel of normals site. In that case, Mutect2 proceeds with assembly and genotyping but marks the panel of normals site for filtering by FilterMutectCalls.\

    here I mean the column of FILTER and INFO, you can see, there are some label have no corresponding value in the INFO column.

    The table you posted below from the Mutect2 documentation answers this question. For example, the MBQ key in the VCF is used for base quality filtering.

    if it is ok. hope we can have a talk in skype, thanks a lot.

    I'm too busy improving Mutect2, but feel free to ask specific questions, one per thread, with the command line used and example output, on the GATK forum.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    The germline resource gatk offered, the_ chromosome name does not start with "chr"_, we need to convert it by myself. the following command is ok

    This belongs on its own thread. Please include the command lines you used for Mutect2 and FilterMutectCalls, along with example output of each.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    so can you help me finding out the relationship between likelihood and count in this?

    @picard_gatk_mj The Mutect2 likelihoods are explained here: https://github.com/broadinstitute/gatk/blob/master/docs/mutect/mutect.pdf. They do presuppose some knowledge of statistics and probability.

  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018

    @davidben said:

    The germline resource gatk offered, the_ chromosome name does not start with "chr"_, we need to convert it by myself. the following command is ok

    This belongs on its own thread. Please include the command lines you used for Mutect2 and FilterMutectCalls, along with example output of each.

    _
    Thanks a lot, Thank you very much for helping me during the holidays._

    creat pon by the function: CreateSomaticPanelOfNormals (so there is no problem in my Pon.vcf)

    here is my command (both have two to compare the results with or without vcf)

    Mutect2
    gatk Mutect2 -L abc.bed -I NGS123_sorted.recal.bam -R hg19.fasta -tumor NGS123 -pon PoN.vcf --native-pair-hmm-threads 10 -O NGS123_pon_filter_somatic.raw.vcf

    gatk Mutect2 -L abc.bed -I NGS123_sorted.recal.bam -R hg19.fasta -tumor NGS123 --native-pair-hmm-threads 10 -O NGS123_somatic.raw.vcf

    FilterMutectCalls
    gatk FilterMutectCalls -V NGS123_pon_filter_somatic.raw.vcf -O NGS123_pon_filter_oncefiltered.vcf

    gatk FilterMutectCalls -V NGS123_somatic.raw.vcf -O NGS123_oncefiltered.vcf

    CollectSequencingArtifactMetrics
    gatk CollectSequencingArtifactMetrics -I NGS123_sorted.recal.bam -R hg19.fasta -O NGS123__tumor_artifact --FILE_EXTENSION \".txt\"

    _FilterByOrientationBias _
    gatk FilterByOrientationBias --artifact-modes G/T -V NGS123_oncefiltered.vcf -P NGS123_tumor_artifact.pre_adapter_detail_metrics.txt -O NGS123_twicefiltered.vcf

    gatk FilterByOrientationBias --artifact-modes G/T -V NGS123_pon_filter_oncefiltered.vcf -P NGS123_tumor_artifact.pre_adapter_detail_metrics.txt -O NGS123_pon_filter_twicefiltered.vcf

    I compare the vcf sites of NGS123_pon_filter_twicefiltered.vcf and NGS123_twicefiltered.vcf.

    I can not undertand your answer_** "Sometimes an assembly region triggered by a site that is not in the panel of normals will include a panel of normals site. In that case, Mutect2 proceeds with assembly and genotyping but marks the panel of normals site for filtering by FilterMutectCalls.\"**_.

    as shlee ever said in https://software.broadinstitute.org/gatk/blog?id=10911.

    Filtering of sites in the panel of normals (PoN) and the matched normal remains unchanged, except that the tool will prefilter most of these such that site records are absent from the VCF.

    so it is hard for me to undertand and distinguish whether the site triggered reassembly(I guess when you said assemble means reassembly) is in pon or not.
    or to be more clear, can you clarify what shell said (filter the sites in pon), for example, give the condition when the sites filtered(you asked what dropped mean, dropped means filter automatically, not labeled with panle_of pon), when they will be labeled.

    The following are the results, I will give a detailed description about it.(I will annotate as F1, F2,F3...)
    F1 : some sites that filter automatically by PON,means they do not exsit in NGS123_pon_filter_twicefiltered.vcf.

    F3: many samples labled with panel_of_pon, not automatically filter,_ so it is contradictory with that "Filtering of sites in the panel of normals (PoN) and the matched normal remains unchanged"_ or filter here means a lot ?

    POP_AF=1.000e-03;P_GERMLINE=-2.597e-04.
    why P_GERMLINE is minus, it is posterior probability as doc said.
    and I do not use argument germline_resoure, why here has POP_AF, and the value is a constant, while if I use germline_resource, _there is no variant site in my vcf after the command _FilterMutectCalls_, I mean ther is not a site at all, _
    and the vcf after Mutetc2 command the POP_AF is what I set ,just like 0.0000025.

    F2: the sites used in F1 and F3 exists in Pon.vcf. also focus on why some labeled, some filter automatically?

    Here is another new sample.
    F4: a new site appear after pon filter(startle),and this site is like follows in the vcf (here vcf means from many normal samples, the eight column vcf, you know my meaning)

    F5 : there is only site in the Pon.vcf of all the 7 vairant sites

    Thanks a lot.

    Post edited by picard_gatk_mj on
  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    I can not understand your answer_** "Sometimes an assembly region triggered by a site that is not in the panel of normals will include a panel of normals site. In that case, Mutect2 proceeds with assembly and genotyping but marks the panel of normals site for filtering by FilterMutectCalls.\"**_.

    Let me give a specific example. Suppose there are two sites in chr17, one at 1234567 that is in the PoN (let's suppose it's an artifact) and another at 1234569 that is not (let's say it's a real somatic mutation). When Mutect2 sees the first site it will realize that it is in the PoN and do nothing, but when it reaches the second site it will see evidence to begin local assembly. Since assembly is done over a region, not just a single site, the assembly region will include the first site. Thus the apparent "variation" of the first site will be assembled as a potential haplotype and genotyped. Mutect2 will mark it with the PON vcf info field, and then FilterMutectCalls will add the panel_of_normals filter.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    Your boxes F1, F2, and F3 are consistent with the documented behavior of Mutect2.

    As for F4 and F5, I need to preface what I'm about to write with the disclaimer that because this is a very filtered, almost certainly artifactual, call, this is a purely pedantic discussion. I'll try to give a general understanding of why this happens, but in the end we should all be aware that it doesn't really matter.

    Anyway, what I think is happening is that you have calls at 9028, 9063, 9072 (using only the last four digits for clarity). Of these, the evidence for 9072 is weak, the evidence for 9028 is very weak, and the evidence for 9063 is strong, but it's in the PoN. Since Mutect2 determines assembly regions in terms of nearby sites with evidence of variation, the effect of the variant at 9063 is to stretch the assembly region containing 9072 to the left. This doesn't happen when you run Mutect2 with the PoN, because 9063 is not considered to exhibit variation.

    Thus: When you use the PoN, site 9028 is not assembled along with 9072 and gets its own separate assembly region. When you don't use the PoN, they are assembled together. Since 9028 is so marginal and very, very near the tumor lod threshold to emit, the fact that in one case it's in the middle of its own assembly region and in one case at the tail of 9072's assembly region is enough to shift the Pair-HMM likelihoods very slightly, probably from 3.3 to 2.9 or something.

    You can test my theory by running Mutect2 with tumor-emission-lod set to something lower, like 1.0. I suspect that you will now find the variant at 9028 in both vcfs, with slightly different TLODs.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    and I do not use argument germline_resoure, why here has POP_AF, and the value is a constant, while if I use germline_resource, there is no variant site in my vcf after the command _FilterMutectCalls, I mean there is not a site at all, _

    In tumor-only mode Mutect2 uses a default population allele frequency. This was too high in 4.0.0.0 and led to too much germline_risk filtering. The solution is to use the latest release of Mutect2.

  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018

    @davidben. thanks a lot.
    I know why you are so busy now, because you know so much, but I know nothing about it, you need a lot of time to acquire further knowledge about it, I feel ashamed of myself.

    but about answer
    "
    This doesn't happen when you run Mutect2 with the PoN, because 9063 is not considered to exhibit variation.

    Thus: When you use the PoN, site_ 9028 i_s not assembled along with 9072 and gets its own separate assembly region. When you don't use the PoN, they are assembled together
    "
    should here 9028 changed to 9063 when you said "site_ 9028 i_s not assembled along with 9072 and". because 9063 appears as a new site.

    "
    Since Mutect2 determines assembly regions in terms of nearby sites with evidence of variation, the effect of the variant at 9063 is to stretch the assembly region containing 9072 to the left. This doesn't happen when you run Mutect2 with the PoN, because 9063 is not considered to exhibit variation.

    the answer is a little hard to understand.
    and you said 9063 has the the biggest signal, if it is not considered to be a variant, 90672 and 9028 certainly not, and I think 9063 maybe assembly in a region, which not containing 9082 and 9072 when run mutect2 with pon, do you think it is maybe a better answer ? but it is a little funny, the same site, it should be the same assembly way with or without pon.
    and when with the PoN, 9063 is not considered to exhibit variation. so it is more hard to believe it should appear in the final vcf, I mean it should filtered automatically. because the left(9028) and right region(9072) of 9063 is weaker,
    "
    and why new site appears, it is very abnormal according to your experience?

    and it is really sad there is no "PASS " in the filter, means there is no good variant.

    Post edited by picard_gatk_mj on
  • picard_gatk_mjpicard_gatk_mj Member
    edited December 2018

    you told me to use the latest, but it has such words when running.
    "
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    Warning: FilterByOrientationBias is an EXPERIMENTAL tool and should not be used for production

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    A USER ERROR has occurred: min-base-quality-score is not a recognized option
    "

    4.0.12.0 seems to not be a production version?
    thanks a lot

  • the POP_AF=1.000e-03 in gatk4.0.0.0, but POP_AF=5.000e-08 in gatk4.0.11.0, really a big change, after use this value, there are some "PASS" in the column of FILTER, what make your team make so big change, whether your team find out some backgroud theory.

    if I use the germline_resource, and set --af-of-alleles-not-in-resource 0.0000025 , so this value will overwrite the value when it is default tumor-only mode , am I right?

    so how did the POP_AF value afftect the Mutect2 apply reassmebly and the final choice of whether this is variant or not.

    thanks a lot

  • another important question is
    you know there are some sites we are very concerned, or whitelist, we want to know they to call even they will be filtered.

    I see there are four argument about this
    "
    --genotype-germline-sites
    false (EXPERIMENTAL) Call all apparent germline site even though they will ultimately be filtered.

    --genotype-pon-sites
    false Call sites in the PoN even though they will ultimately be filtered.

    --genotyping-mode
    DISCOVERY Specifies how to determine the alternate alleles to use for genotyping

    --alleles
    null The set of alleles at which to genotype when --genotyping_mode is GENOTYPE_GIVEN_ALLELES
    "

    I guess --alleles and --genotyping-mode is a pair for use according to the literal surface,
    but how they I give the value of --alleles, can you give me an example. it will save a lot of time if we just cared about some gene or some sites of gene, so can the value set to gene level or set leve

    when set --genotype-pon-sites and --genotype-germline-sites, if we put aside the time cost, is there any accuracy or site numbers difference?

    thanks a lot.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    @picard_gatk_mj

    should here 9028 changed to 9063 when you said "site_ 9028 i_s not assembled along with 9072 and". because 9063 appears as a new site.

    In your screenshot above 9028 is the new site. The answer above does not need to be changed.

    and I think 9063 maybe assembly in a region, which not containing 9082 and 9072 when run mutect2 with pon, do you think it is maybe a better answer ? but it is a little funny, the same site, it should be the same assembly way with or without pon.

    I think the answer I gave above is the correct one, but as I mentioned above you can test it by lowering the emission lod threshold.

    and why new site appears, it is very abnormal according to your experience?

    It's uncommon enough that I had to think for a while to come up with the explanation i gave above.

    and it is really sad there is no "PASS " in the filter, means there is no good variant.

    Perhaps it is sad, but sequencing artifacts are common and real somatic mutations are not.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    you told me to use the latest, but it has such words when running.

    Those words were present in 4.0.0.0, but FilterByOrientationBias is a perfectly reliable tool, although we now recommend the new orientation bias tool.

    A USER ERROR has occurred: min-base-quality-score is not a recognized option

    --min-base-quality-score is an argument for Mutect2, but it is not an argument for FilterByOrientationBias. This error message is to be expected.

    4.0.12.0 seems to not be a production version?

    It's a production version, as is every release.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    the POP_AF=1.000e-03 in gatk4.0.0.0, but POP_AF=5.000e-08 in gatk4.0.11.0, really a big change, after use this value, there are some "PASS" in the column of FILTER, what make your team make so big change, whether your team find out some backgroud theory.

    This is described in the Mutect2 documentation: https://github.com/broadinstitute/gatk/blob/master/docs/mutect/mutect.pdf

    if I use the germline_resource, and set --af-of-alleles-not-in-resource 0.0000025 , so this value will overwrite the value when it is default tumor-only mode , am I right?

    Specifying arguments on the command line overrides default values, yes.

    so how did the POP_AF value afftect the Mutect2 apply reassmebly and the final choice of whether this is variant or not.

    This is also described in the Mutect2 documentation. The population allele frequency is used in the germline model and the contamination model.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    I guess --alleles and --genotyping-mode is a pair for use according to the literal surface,

    but how they I give the value of --alleles, can you give me an example. it will save a lot of time if we just cared about some gene or some sites of gene, so can the value set to gene level or set level.

    The argument for --alleles is a vcf of alleles to forcibly genotype, regardless of how little evidence there is. Ordinary filtering will still occur, but you are guaranteed to see these sites in the filtered Mutect2 / FilterMutectCalls output. If you want to override filtering you can always do something like

    gatk SelectVariants -V filtered_mutect_calls.vcf --concordance whitelist.vcf -O whitelisted_mutect_calls.vcf
    

    when set --genotype-pon-sites and --genotype-germline-sites, if we put aside the time cost, is there any accuracy or site numbers difference?

    Aside from the germline and pon sites that would otherwise not be genotyped, no.

  • @davidben, thanks a lot. do you know why the P_GERMLINE is minus?

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    It used to be log scale, but now it's phred-scaled.

  • log10(Posterior probability ) (used to )
    changed to
    -10*log10(Posterior probability ) (now)

    used to means gatk version until which.

    thanks, @davidben

Sign In or Register to comment.