Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

tumor-only mode, should CollectSequencingArtifactMetrics and FilterByOrientationBias be used on bam?

I use the latest version of gatk4 and when in tumor-only mode, should I do CollectSequencingArtifactMetrics and FilterByOrientationBias be used on sorted.recal.bam ?
because in paired mode, CollectSequencingArtifactMetrics and FilterByOrientationBias need to be only used in normal bam,

so if need to on tumor bam why? if not, why, thanks a lot

here I also want to ask another question, when make pon of snv call, should the following three arguments be used?
--germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz \
--af-of-alleles-not-in-resource 0.0000025 \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \

be used to make the each sample vcf? and if so, should the value of --af-of-alleles-not-in-resource change a little?
I see in your tutorial (How to) Call somatic mutations using GATK4 Mutect2 , you just used --disable-read-filter to generate vcf for pon, can you clarify here

gatk Mutect2 \
-R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta \
-I HG00190.bam \
-tumor HG00190 \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
-L chr17plus.interval_list \
-O 3_HG00190.vcf.gz



and we really hope gatk staff can liftover a _af-only-gnomad_hg19.vcf.gz for us, because we really are afraid making it wrong.

by the way, will --germline-resource make a big contribution for germline variant filtering

I am also a little confused about the following sentence in your tutorial (How to) Call somatic mutations using GATK4 Mutect2

If a variant is absent from a given germline resource, then the value for --af_of_alleles_not_in_resource applies. For example, gnomAD's 16,000 samples (~32,000 homologs per locus) becomes a probability of one in 32,000 or less. Thus, an allele's absence from the germline resource becomes evidence that it is not a germline variant.

here it seems to just rely on if the variant is in the a given germline resource, then discard the variant, and I do not see any significance of --af-of-alleles-not-in-resource 0.0000025 here ,

Thanks a lot, thanks for all the staff in gatk , they helped me so much, wish they all a better future.

Answers

  • 29043594952904359495 Member
    anyone helps, thanks a lot
  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    @2904359495

    because in paired mode, CollectSequencingArtifactMetrics and FilterByOrientationBias need to be only used in normal bam

    It's the other way around -- it only needs to be run on the tumor sample(s), and you should run it any time you suspect sample quality issues such as with FFPE samples.

    here I also want to ask another question, when make pon of snv call, should the following three arguments be used?
    --germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz \
    --af-of-alleles-not-in-resource 0.0000025 \
    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter

    You should not use any of these. --af-of-alleles-not-in-resource 0.0000025 should now only be used if you have a really good reason to override the default. --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter is now the default. As for the germline resource, don't use it in Mutect2 when making a pon, but you may use it as an argument to CreateSomaticPanelOfNormals as of the 4.1.1 release.

    and we really hope gatk staff can liftover a _af-only-gnomad_hg19.vcf.gz for us, because we really are afraid making it wrong.

    It's in the public GATK resource bucket: gs://gatk-best-practices/somatic-hg38/af-only-gnomad.hg38.vcf.gz

    by the way, will --germline-resource make a big contribution for germline variant filtering

    Yes, a huge difference. Don't run Mutect2 without gnomAD, except for making a pon. And mouse sample etc.

    here it seems to just rely on if the variant is in the a given germline resource, then discard the variant

    Mutect2 does not simply hard filter based on preence in the germline resource. It uses the population allele frequency in a probabilistic model.

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    One more thing: the old orientation bias tool is now deprecated. We encourage you to use LearnReadOrientationModel now. A blog post will come out within a coupe of weeks.

  • 29043594952904359495 Member
    @davidben
    Thanks a lot. I am really deeply touched when you answer me,

    First I want to confrim somethings with you.
    No matter in tumot-only or matched sample, CollectSequencingArtifactMetrics only need to do on tumor sample.

    secondly, if I use argument --germline-resource, there is not a must to set --af-of-alleles-not-in-resource value, because it has a default value if I use argument --germline-resource

    thirdly, It's in the public GATK resource bucket: gs://gatk-best-practices/somatic-hg38/af-only-gnomad.hg38.vcf.gz, but there is no hg19, I hope you can supply that, through crossmap I may can do, but I am afraid it is wrong .



    I have read the pdf of gatk wokkshop of 201903, the LearnReadOrientationModel command result now passed to Mutect2 .

    gatk CollectF1R2Counts \
    -R GRCh38.fasta \
    -I tumor.bam \
    -O tumor-artifact-prior-table.tsv \
    -alt-table tumor-alt.tsv \
    -ref-hist tumor-ref.metrics \
    -alt-hist tumor-alt.metrics

    gatk LearnReadOrientationModel \
    -alt-table my-tumor-sample-alt.tsv \
    -ref-hist my-tumor-sample-ref.metrics \
    -alt-hist my-tumor-sample-alt-depth1.metrics \
    -O my-tumor-sample-artifact-prior.tsv

    gatk Mutect2 \
    -R ref_fasta.fa \
    -I sample.bam \
    -L intervals.list \
    -O sample.vcf.gz \
    --orientation-bias-artifact-priors tumor-sample-artifact-prior.tsv







    But there is no a clear documentation of how the argument value made, I tried as the pdf, but it has a error, I used gatk4.1.1.0,
    /public/software/gatk-4.1.1.0/gatk CollectF1R2Counts -R hg19/ucsc.hg19.fasta -I NGS190124-10DE.sorted.recal.bam -O tumor-artifact-prior-table.tsv -alt-table tumor-alt.tsv -ref-hist tumor-ref.metrics -alt-hist tumor-alt.metrics,

    but it get an error "A USER ERROR has occurred: a is not a recognized option

    ***********************************************************************
    Set the system property GATK_STACKTRACE_ON_USER_EXCEPTION (--java-options '-DGATK_STACKTRACE_ON_USER_EXCEPTION=true') to print the stack trace."

    I also do not find argument -alt-table and -ref-hist and -alt-hist in CollectF1R2Counts.

    here I supply a figure,
  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    No matter in tumot-only or matched sample, CollectSequencingArtifactMetrics only need to do on tumor sample.

    Yes, that's right.

    secondly, if I use argument --germline-resource, there is not a must to set --af-of-alleles-not-in-resource value, because it has a default value if I use argument --germline-resource.

    It always has a default value. --af-of-alleles-not-in-resource is only for users with a very, very good reason not to use the default.

    thirdly, It's in the public GATK resource bucket: gs://gatk-best-practices/somatic-hg38/af-only-gnomad.hg38.vcf.gz, but there is no hg19, I hope you can supply that, through crossmap I may can do, but I am afraid it is wrong .

    There is a somatic-b37 folder, which is equivalent to hg19.

    I tried as the pdf, but it has a error, I used gatk4.1.1.0

    We optimized this part of the workflow in the 4.1.1 release and changed the parameters. A blog post is coming out soon, but in the meantime here's a draft:

    There are three steps to the filter. First, run Mutect2 with the --f1r2-tar-gz argument. This creates an output with raw data used to learn the orientation bias model. This works in multi-sample mode as well. That is, when multiple tumor samples are specified, you only need a single --f1r2-tar-gz output, which contains data for each tumor sample.

    gatk Mutect2 -R ref.fasta \
       -L intervals.interval_list \
       -I tumor.bam \
       [-germline-resource af-only-gnomad.vcf] \
       [-pon panel_of_normals.vcf ]  \
       --f1r2-tar-gz f1r2.tar.gz \
       -O unfiltered.vcf
    

    Next, pass this raw data to LearnReadOrientationModel:

    gatk LearnReadOrientationModel -I f1r2.tar.gz -O read-orientation-model.tar.gz
    

    Finally, pass the learned read orientation model to FilterMutectCallswith the -ob-priors argument:

    gatk FilterMutectCalls -V unfiltered.vcf \
       [--tumor-segmentation segments.table] \
       [--contamination-table contamination1.table] \
       --ob-priors read-orientation-model.tar.gz \
       -O filtered.vcf
    
  • 29043594952904359495 Member

    @davidben Thanks a lot for your continous help, I think gatk really do a lot to improve the accuracy, we really appreciate your big contribution.

    here I still want to confirm somethings with you

    1 CollectF1R2Counts has been integrated into mutect2 , and there is separate CollectF1R2Counts command or at least there no need to do separate CollectF1R2Counts .

    2 you also mentioned a huge message "This works in multi-sample mode as well.", so do you mean mutect2 now support many tumor samples and one matched normal sample, just a little similar to genetype in germline calling, this is really huge message, because we can usually get Primary focus and Metastatic Tumor of a same patient, but just one matched blood, so is there a example in gatk documentation?

    3 here is really a big question, because many people have asked geomAD for hg19, nearly all the staff in gatk told them to do liftover, but no people ever said here is a somatic-b37 folder, which is equivalent to hg19.,
    I looked the folder, the files size differs greatly with hg38, I guess it is because it is raw though I do not know what is raw means

    thanks a lot, you are really a very good person

  • 29043594952904359495 Member

    another question is
    gatk GenomicsDBImport -R reference.fasta -L intervals.interval_list \
    --genomicsdb-workspace-path pon_db \
    -V normal1.vcf.gz \
    -V normal2.vcf.gz \
    -V normal3.vcf.gz

    gatk CreateSomaticPanelOfNormals -R reference.fasta -V gendb://pon_db -O pon.vcf.gz

    what gendb://pon_db means , it this a file, or it is a special database, and whether4.1.1.0 still support
    gatk CreateSomaticPanelOfNormals \
    -vcfs 3_HG00190.vcf.gz \
    -vcfs 4_NA19771.vcf.gz \
    -vcfs 5_HG02759.vcf.gz \
    -O 6_threesamplepon.vcf.gz

    is the new command more accurate?

  • 29043594952904359495 Member

    @davidben ,thanks a lot, please help me when you are not thta busy

  • 29043594952904359495 Member
    edited April 8

    @davidben I think liftover may not be a good choice, because variant in one hg38 genome version maybe not in hg19, so it is better to use the origin hg19, so it is convenient for you to supply this.

    what is more, in the gnomAD, it clearly said it only has hg19
    What genome build is the gnomAD data based on?
    All data are based on GRCh37/hg19.

    you can see the geomAD just supply hg19, so you must have get hg19 first, then you liftover it to hg38

  • 29043594952904359495 Member
    edited April 8

    @davidben
    As for the germline resource, don't use it in Mutect2 when making a pon, but you may use it as an argument to CreateSomaticPanelOfNormals as of the 4.1.1 release.

    so do you recommend use germline resource for CreateSomaticPanelOfNormals command, I am afraid it may infect so much, so I want to ask for your advice.

    another question is still about creat pon, it the 201903 workshop, the result of LearnReadOrientationModel passed to Mutect2 and the result of CalculateContamination passed to FilterMutectCalls.

    and according to your latest answer, both result now passed to FilterMutectCalls, am I right?

    so when make normal sample vcf, there is no addtional parameters to add for Mutect2?

    just like so plain command

    gatk Mutect2 \
    -R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta \
    -I HG00190.bam \
    -tumor HG00190 \
    -L chr17plus.interval_list \
    -O 3_HG00190.vcf.gz
    
    gatk CreateSomaticPanelOfNormals \
    -vcfs 3_HG00190.vcf.gz \
    -vcfs 4_NA19771.vcf.gz \
    -vcfs 5_HG02759.vcf.gz \
    -O 6_threesamplepon.vcf.gz
    

    and when do sample analysis, add the LearnReadOrientationModel and CalculateContamination command, and pass both the file to FilterMutectCalls, am I right?

    by the way, is the newest blog about how to use 4.1.1.0 generated?

    thanks a lot

    Post edited by 2904359495 on
  • 29043594952904359495 Member

    @SkyWarrior do you have some suggestions, thanks a lot

  • 29043594952904359495 Member
    edited April 9

    @davidben
    gatk CreateSomaticPanelOfNormals \
    -vcfs 3_HG00190.vcf.gz \
    -vcfs 4_NA19771.vcf.gz \
    -vcfs 5_HG02759.vcf.gz \
    -O 6_threesamplepon.vcf.gz

    in 4.1.1.0 , this is not a write code , it will raise error "A USER ERROR has occurred: v is not a recognized option, it mistake -vcf as -v, so we must use like the newest code?

    Step 1. Run Mutect2 in tumor-only mode for each normal sample.
    gatk Mutect2 -R reference.fasta -I normal1.bam -O normal1.vcf.gz

    Step 2. Create a GenomicsDB from the normal Mutect2 calls.
    gatk GenomicsDBImport -R reference.fasta -L intervals.interval_list \
    --genomicsdb-workspace-path pon_db \
    -V normal1.vcf.gz \
    -V normal2.vcf.gz \
    -V normal3.vcf.gz

    Step 3. Combine the normal calls using CreateSomaticPanelOfNormals.
    gatk CreateSomaticPanelOfNormals -R reference.fasta -V gendb://pon_db -O pon.vcf.gz

    so I run the step 2,
    /gatk GenomicsDBImport -R ucsc.hg19.fasta -L TT.bed \
    --genomicsdb-workspace-path /xx/xx/pon_db \ /ABCA06.somatic.raw.vcf -V .......

    it raise another error
    A USER ERROR has occurred: Bad input: GenomicsDBImport does not support GVCFs with MNPs. MNP found at chr1:11182061 in VCF H02.somatic.raw.vcf

    and the site information looks like following(really a mnp)
    chr1 11182061 . TT GC . . DP=2;ECNT=1;MBQ=0,10;MFRL=0,0;MMQ=60,60;MPOS=27;POPAF=7.30;TLOD=6.34 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:0,2:0.750:2:0,0:0,0:0,0,0,2

    I think mnp is a very common question, you must have thougth that question,
    but the pon_db still gererated.
    so I tried to run step3
    gatk CreateSomaticPanelOfNormals -R reference.fasta -V gendb:/xx/xx/pon_db -O /pon.vcf.gz

    but still raise another error

    so help me please, i doubt is this command is still in beat, can not put into production, and is 4.1.1.0 a version can be put into production ?
    thanks a lot

  • syer89syer89 Member
    edited May 20
    Hi @davidben thanks for laying out the three steps to filter for the newest version as there is no clear documentation anywhere yet.
    Does this mean FilterByOrientationBias (and CollectSequencingArtifactMetrics) need not be done after FilterMutectCalls (Tutorial#11136) ?
    Also, I haven't seems any clear doc regarding FilterAlignmentAritifacts. Any recommendations when it needs to be used and before/after FilterMutect calls ? Thanks for all your help!
  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    @syer The new Mutect2 tutorial describing the new read orientation bias model should be coming out any day now. Please stay posted for it. CollectSequencingArtifactMetrics and FilterByOrientationBias are now deprecated and should not be used. FilterAlignmentAritifacts is not documented because it's not part of our best practices. It sacrifices too much sensitivity. We are working on an improved version that realigns an assembled unitig of alt-supporting read pairs rather than realigning each pair in isolation. It's an appealing idea and I hope it works but we can't make any promises yet.

  • syer89syer89 Member
    Thanks @davidben. So to confirm, >=GATK v4.1.1 for now we just stick with the 3 steps you mentioned for filtering with mutect2.
  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    That's correct.

Sign In or Register to comment.