The DP value of variants called in Mutect2 vcf much lower than depth in input bam file
a) GATK version used: 4.6.0.0
Input: Custom targeted panel sequencing, hybridization capture based, brain tissue samples, average deduplicated sequence depth ~1000-1500X
I am using Mutect2 in GATK 4.6.0.0 to call indels and SNVs. We have done all proper pre-processing (fastqc, alignment to ref genome with bowtie2, removing duplicates with picard). The vendor who sold us the library prep kit confirmed that the input sequencing data is of good quality with a >70% on-target rate. The vendor who did the sequencing confirms that sequencing went well. I am therefore confused as to why we start with a bam with average depth of ~1300X, but the output mutect2 file only has an average depth (DP) of ~100-300X.
In reading previous similar forums, I wonder if maybe downsampling could be contributing to this, but I read that that usually applies to amplicon-based sequencing, and we used hybridization capture. Are there other reasons why the depth for called variants in the vcf is so low? I'm new to this kind of analysis, so any assistance would be much appreciated. Thanks!
-
There will be a few questions to ask before we proceed.
1- How did you calculate the average on target depth for these samples? Did you use CollectHsMetrics?
2- Have you tried with --max-reads-per-alignment-start 0 option with Mutect2 ? Mutect2 and HaplotypeCaller uses downsampling to ease the performance however for those who are absolutely in need for getting every single read to be processed we recommend using this option regardless of the data type.
Finally we do not recommend using bowtie2 for aligner with Mutect2 unless you have a very important reason to do so. Bowtie2 uses a different scoring scheme for MAPQ which is not compatible with GetPileupSummaries tool that we have as an accessory for somatic calling.
Let us know the answers to these questions and we can proceed further.
Regards.
-
Hi Gökalp Çelik!
Thanks for the quick reply and your advice.
1. For calculating the depth, we used the depth operator in samtools with a path to our covered .bed file and the .bam file after bowtie2 and then the same with the .bam file after mutect2
2. I also feared downsampling and set the downsampling stride to 20 in a recent run of a few samples. It seems like depth grossly improved. Is this similar to what you propose? Should I employ this max-reads-per-alignment method instead?
3. We started using bowtie2 at the recommendation of another lab who used it for their pipeline. I think their only rationale was that it was already installed on our computing cluster. If not bowtie2, what would you recommend? BW2? Some other alignment tool?
Thanks again for all your wisdom!
-
Hi again.
You can use --max-reads-per-alignment-start 0 parameter to disable downsampling in Mutect2.
Bowtie2 uses maximum mapping quality of 44 for local mode and 42 for end-to-end mapping. These values are not directly compatible with GetPileupSummaries tool which uses minimum MAPQ 50 therefore if you use the tool without modifying this parameter it will return an empty file. If this tool is not of your interest then it may be OK to proceed. Overall BWA MEM can map more reads compared to Bowtie2 and may recover more shorter matches.
I hope this helps.
Regards.
-
Gökalp Çelik Hey there again! Thank you again for this helpful information. One last question just to confirm: turning off downsampling doesn't have any deleterious implications for downstream analysis, right? The reads will still be high quality, not duplicates, etc?
-
No it does not affect the way other ReadFilters perform. It will only disable subsampling of reads.
Regards.
Please sign in to leave a comment.
5 comments