Holiday Notice:
The Frontline Support team will be offline February 18 for President's Day but will be back February 19th. Thank you for your patience as we get to all of your questions!

which file to use as input to DepthOfCoverage?

mayaabmayaab IsraelMember ✭✭

Hello GATK team,

I'm running the Best Practice pipeline, and I would like to do some coverage calculation on my samples.
I want to use DepthOfCoverage, but I can't decide which bam file to use - the original one, the one after BQSR stage (right before HC) or the bamout of HC?
The point is that the first two files don't contain the reads as they were mapped in the variant calling, since the HC changes the mapping dramatically. The bamout contains haplotypes, and therefor reads are merged together, and the coverage might look lower than the one that uses the reads.

What do you recommend?
Maya

Tagged:

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Maya,

    It depends what you want to do with the coverage information. If you just want the raw initial coverage for general QC, then use the original bam file. If you want to know what is the coverage that is actually usable for variant calling, you can use the bamout file. It contains all the reads that were deemed informative and were effectively used by HaplotypeCaller. Any other reads are useless, so including them in your coverage counts makes it look like you have more data than you can actually use.

    Better yet, do both and run a comparison, which will tell you, out of the raw sequence that you paid for, how much was actually usable. If you don't like the resulting ratio, complain to your sequencing provider.

  • mayaabmayaab IsraelMember ✭✭

    The original bam contains ~41 million reads, and the bamout contains ~4 million reads.
    I don't think that the number of usable reads is ~4 million. Does the bamout indeed contains the informative reads? doesn't it contains the haplotypes created by the informative reads?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @mayaab Did you read the documentation about what -bamout does? And did you try opening the file and looking at it in a genome browser? Those are the first things you should think of doing in a situation like this.

    The file contains the haplotypes and the informative reads for the active regions that were processes. It does not contain any data for regions that were considered inactive (non variant).

  • mayaabmayaab IsraelMember ✭✭

    I understand this. I'm asking this question because you said:
    "Better yet, do both and run a comparison, which will tell you, out of the raw sequence that you paid for, how much was actually usable"
    but if the bamout contains the haplotypes. AS far as I understand the coverage calculated on the haplotypes is not comparable to the coverage of the original bam, since the haplotypes marge overlap reads.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    What I'm trying to suggest is that you can open up the bamout in IGV and do a sanity check of how many reads you see there. There will be only two artificial haplotype reads and they will be colored differently than the rest by default, so you should be able to see immediately if there are just a few real reads at the site of interest, or hundreds.

Sign In or Register to comment.