Hello, all! I'm trying to understand the nuances between the MosDepth Summary Statistic 'Bases' and the Samtools Summary Statistic 'Mapped Bases' in the General Stats. I'm getting pretty vastly different values from the two tools as seen in the scatterplot attached (note the scales are different orders of magnitude) What I'm trying to do is get information on how much of the genome was actually covered in each sample. I ran the pipeline without a regions file as this is an amplicon experiment with degenerate primers, so we don't know the extent of what they'll cover until we dig into the data. Therefore, I just ran everything as if it were WGS. Does the fact that MosDepth does 500 bp windows have something to do with this? How would you suggest I properly / fairly assess how much of the genome we've covered? Alternatively, I can create a bed file from the .cram alignments and actually see how many bases are covered, too, if neither of the metrics I've mentioned are what I want. Thanks for any advice 🙂

It's a bit strange that the Samtools mapped bases differ so much from Mosdepth total bases, but they are expected to differ quite a bit by the way alignments are filtered, if mate overlaps and read clippings are counted, etc. Total bases is independent of the window size in Mosdepth, which only affects the regions for which coverages are reported. If you want to know how much of the genome is covered at a certain depth threshold, you can check the

*.mosdepth.global.dist.txt

in the Mosdepth output.

In addition to the mate overlaps, etc., Alexander mentioned, I wonder how multimapping reads are being treated.

samtools

tends to be take the simplest possible approach, e.g.

samtools view -c

counts basically the lines in the SAM/BAM -- it does not care about read pairs or multimappers, so if you wanted to count fragments, you would have to apply various filters or other tricks to avoid double-counting. I would expect MosDepth to be more sophisticated and "intelligent".

Thank you both for your input! Unfortunately, the

*.mosdepth.global.dist.txt

isn't helpful here because I don't have a defined regions file, so it is passing in whole chromosomes and reporting 0x coverage for each. The more I think about this and how MosDepth works, I think it's reflecting the amount of unique/nonduplicated genome coverage per sample. I did some investigating to see what other correlations I could find in the general stats, and found that if I plot Unique Reads (Total Reads - Duplicates) vs SNPs discovered, I see a positive correlation, similar to what I saw with the MosDepth coverage (scatterplot attached). I think I will try to investigate the output .cram files from markduplicates to see how much genomic space is covered. Do you think that's a good idea?

Mosdepth by default filters any reads with any of the bits in the 1796 flag set (unmapped, not primary, failed QC, PCR or optical duplicate). Samtools stats seems to only care about the mapped bit. A possible reason for the difference could therefore be the secondary mappers and PCR duplicates that are counted by Samtools.

I likely need to redo some settings as my global coverage file only reports on 0x, 1x, or 2x coverage, so everything's just getting binned into 0x. Here's what the file looks like for one sample

I think your problem is that only a tiny fraction of your genome is covered, so the coverage at depth>=1 is rounded to 0.00. A better approach in this case is to look at depth either in small windows or with basepair-resolution, depending on the size of your reference.

Yeah, I definitely agree. We weren't sure before getting these sequencing results what fraction of the genome would be covered at 1x which is why I'm trying to explore which parts of the genome we actually covered. We knew it would be a small fraction, but we didn't know how small. Based on our discussion I think the MosDepth coverage is a better surface level reflection of true coverage as compared to Samtools. Now I just need to use other tools or refine my strategy (with smaller bins) to get a real number for what fraction of the genome we captured. Hopefully that clarifies why I'm asking what I'm asking and my takeaways from our conversation make sense.

I think your best bet to get a well defined and useful metric is to calculate it yourself. The

*.per-base.bed.gz

files have the exact coverage (as determined by

mosdepth

with its filtering criteria), and you can intersect this with any regions of interest (e.g. WES targets), if you want, and you can also define a coverage threshold (1x, 5x, 10x, ...) and calculate the fraction of the genome covered at this level. If there is another file that already has numbers that you think are the metric you want, you may still want to verify that from the raw data in

*.per-base.bed.gz

.