Mean and median coverage stats vary heavily between Qualimap and Picard
Here is the difference found for 19 UDN samples.
Activity
Picard's definition:
MEAN_COVERAGE The mean coverage in bases of the genome territory, after all filters are applied. MEDIAN_COVERAGE The median coverage in bases of the genome territory, after all filters are applied.I haven't found qualimap's definition of these metrics. My current guess is that "filters" applied are causing the difference, but not sure why or how.
Which one should we be using in our analyses? I'm leaning picard's but need to have proper reasoning on this.
Edited by Manavalan Gajapathy-
Did some googl'ing but no luck so far. Here is the closest question on this topic but this post doesn't have answer posted - https://www.biostars.org/p/380723/
mentioned in merge request !2 (merged)
I think your correct on why Picard and Qualimap coverage metrics differ. According to Qualimap's docs they note that the meaning of
Xin coverage metrics of10Xis| [1] | Example for the meaning of X: If one genomic region has a coverage of 10X, it means that, on average, 10 different reads are mapped to each nucleotide of the region.coming from the BAM QC section of their docs. This is just a guess but combining that info with the note about the
skip duplicatesoption that statesSkip duplicates This option allows to skip duplicate alignments from analysis. There are three modes of this option. By default, the duplicates are skipped only if they are flagged in BAM file and remaining alignments are further analyzed by Qualimap. Additionally it is possible to skip only the duplicates detected by Qualimap method (based on duplication rate estimation) or aplly both approaches. Number of skipped duplicates will be shown in the report.means that the "filtering" done by qualimap versus Picard is the differentiating factor and that Qualimap's output is closer to "raw" coverage than Picard's is.
With all that said what I'd want to know when answering this question is if the difference between the two outputs would help with identifying samples with potential issues? As in, if we use Picard stats will this highlight samples that will potentially have issues at the variant calling step that otherwise wouldn't be hinted at using the stats from Qualimap?
If I had to guess I would bet the answer is
nobut I have zero tangible data to back that up. If we roll with that assumption though then I'd say we just present the Qualimap stats. The reason being if we have a WGS sample from Vandy that say mean coverage is 30X versus 24x we might be highly inclined to request top-off sequencing for the sample to hit 30X . Again, this operates on the assumption that there is no appreciable difference in variant calls produced. Throw into the mix the differences of what GATK reports as the read depth at a site versus what someone would see in the BAM via IGV and it just becomes easier to explain that Qualimap is the most analogous to "stack the reads, count the bases at that position, average across genome" will make sense to most folks without having to say "there's some special filtering done that isn't quick to explain and that's why we use Picard's filtered metric".-
@bwilk777 I agree with you. I have since moved on to use qualimap's mean coverage, especially when I performed QC review on many projects (post this MR submission). My reasoning was similar to yours as well.
if we use Picard stats will this highlight samples that will potentially have issues at the variant calling step that otherwise wouldn't be hinted at using the stats from Qualimap?
I differ in this though as I think it would highlight bad samples as Picard would filter out mapped reads for various reasons (besides duplicates filtering). However I don't have any data to back this claim.
I am going back and forth on adding Picard's mean-cov among the metrics we check, but it will be more of "keep an eye on this metric" instead of "use it to qualify a sample as pass or fail".