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.

Errors in Mutect2 AD, F1R2 and F2R1 counts

Hello GATK team,

We are using GATK 4.0.10.1, and we noticed error in mutect2 ALT/REF count, that we cannot explain. Some variant will be reported with coverage from forward and reverse read were there is no reverse read with alternative.
We looked at original bam and the --bam-output results with samtools tview and none of them display the same numbers as Mutect2 vcf.

Here is an example with Mutect2 VCF line, and corresponding count in orignal bam and bam-out from mutect2. Mutect2 report 5 alternative in reverse strand where at most we see one in mutect2 bam.

We noticed this for ~2% of driver mutations that where "wrongly(?)" called in our project.

Hopefully it comes from our understanding of the data, if so could you explain to us how to correctly check mutect2 results.

Thanks,

Ismael

############
## Mutect2 VCF
X 76939699 . T C . . DP=27;ECNT=1;POP_AF=5.000e-08;TLOD=32.34 GT:AD:AF:DP:F1R2:F2R1:MBQ:MFRL:MMQ:MPOS:SA_MAP_AF:SA_POST_PROB 0/1:16,9:0.370:25:10,4:6,5:34,42:195,232:60:19:0.222,0.364,0.360:0.069,0.013,0.918

ref = 16
alt = 9

ALT_F1R2: 4
ALT_F2R1: 5
REF_F1R2: 10
REF_F2R1: 6

############
## Original bam

ref = 17
alt = 8

ALT_F1R2: 8
ALT_F2R1: 0
REF_F1R2: 12
REF_F2R1: 5


############
## Mutect2 --bam-output
ref = 18
alt = 9

ALT_F1R2: 8
ALT_F2R1: 1
REF_F1R2: 12
REF_F2R1: 6

Best Answer

Answers

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin
    edited March 19

    Hi @ismael_P

    Please try this with the latest mutect2 GATK4.1 version and report back to us if the results change.

  • ismael_Pismael_P Member
    edited March 20
    Hello,

    I got the same result with GATK 4.1:

    X 76939699 . T C . . DP=27;ECNT=1;MBQ=34,42;MFRL=195,232;MMQ=60,60;MPOS=19;POPAF=7.30;SAAF=0.222,0.364,0.360;SAPP=0.069,0.013,0.918;TLOD=32.3GT:AD:AF:DP:F1R2:F2R1 0/1:16,9:0.370:25:10,4:6,5

    Here is the samtools tview of this position (I'm sorry, I'm not sure it will read well in the forum).
    I count only one "c" on reverse strand (coma), and 8 on forward (dots). So the number of alternative reported by Mutect2 correspond but not the orrientation count.

    76939701 76939711
    TTGGGCACAATTAGTGCGGAA
    Y....................
    ------------------
    C............
    C....................
    .....................
    .....................
    C....................
    .....................
    .....................
    C........ ,,,,,,,,,,,
    .....................
    C....................
    C....................
    ...... ,,,,,,,,,,,,,,
    .....................
    c,,,,,,,,,,,,,,,,,,,,
    ,,,,,,,,,,,,,,,,,,,,,
    ,,,,,,,,,,,,,,,,,,,,,
    ,,,,,,,,,,,,,,,,
    ,,,,,,,,,,,,,,,,,,,,,
    C....................
    .....................
    .....................
    ,,,,,,,,,,,,,,,,,,,,,
    .....................
    C....................
    ,,,,,,,,,,,,,,,,,,,,,
    .....................
    .....................
  • ismael_Pismael_P Member
    Hello,

    I think this is more readable like this. This is the mpileup results from samtools at the variant position.

    Pileup original BAM
    X 76939699 T 24 .CC,CC..,.C.,.C.C..,.C^].^]. CKKDKMCiHCECGCL=KBCFDKE>

    Pileup Mutect2 bam
    X 76939699 T 27 CC..C..C.CC..c,,,,C..,.C,^].^]. KECCLCCKCKMC=BBHDGKBCFDKGE>
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    HI @ismael_P

    Sorry about the delay, I am looking into this now.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @ismael_P

    I have the asked the mutect2 developer to look into this. I will get back to you shortly.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @ismael_P

    Can you please provide IGV screenshots, sorted by base so we see the variant reads, with view as pairs and color by read strand, for the bamout and the original bam?

  • ismael_Pismael_P Member
    Dear @bhanuGandham,

    Here is the IGV views of one problematic position (the same as before) colored by read strand.

    By curiosity, I colored the bamout by "first in pair" and I fall back on the F1R2 and F2R1 counts given by mutect2. In a way this correspond to the F1R2 and F2R1 description : "Count of reads in F1R2 pair orientation supporting each allele".

    If this is the case, I think F1R2 and F2R1 definition are a little bit miss leading and it would be interesting to get the strand coverage of variants in mutect2 output.

    Thank you
  • xiuczxiucz Member

    @davidben , what is the SB VCF field, this didn't occur in gatk4.1.0.0?

  • davidbendavidben BostonMember, Broadie, Dev ✭✭✭

    @xiucz This became a standard Mutect2 annotation in 4.1.1, I think, when we overhauled Mutect2 filtering. An example of the SB annotation is: 23,30,33,18 means that the reference allele is supported by 23 forward reads and 30 reverse reads and the alternate allele is supported by 33 forward reads and 18 reverse reads.

  • xiuczxiucz Member
    edited June 13

    @davidben thank you, I will have a try

Sign In or Register to comment.