Base-by-Base Gene Coverage Analysis in .bam Files Using samtools

[ynkszlr]

Hi everyone!

I’m currently analyzing NGS data and am looking for advice on the best approach to calculate gene coverage on a base-by-base level using .bam files.

Initially, I merged my .bam files and used samtools depth to assess coverage. Then, I mapped the coverage data to specific genes using a .bed file with the relevant coordinates. However, I'm now uncertain if merging all samples like this is a reliable approach, although they were all aligned to the same reference.

To address this, I also tried calculating coverage separately for each sample, averaging the base coverage across samples for each gene. This, however, gave unexpectedly lower and quite different values compared to the merged approach.

My main goal is to assess the gene performance/coverage as a whole, rather than focusing solely on individual sample coverage. I need the base-by-base details to highlight any regions that may lack adequate coverage. Given this, I would like to know:

1. Is merging .bam files and using samtools depth a valid approach for obtaining gene-level coverage?
2. Does averaging base coverage across samples accurately reflect gene coverage, and what factors could cause discrepancies between these two methods?
3. Any recommended best practices for calculating detailed, base-level coverage for specific genes?

Any guidance or shared experience would be really helpful. Thanks in advance!

[fchatonnet]

Hello ynkslzr,

some comments about your question: coverage is for me different from read depth in the way that it should be normalized by total read depth (i.e. library size) for each sample. Having 10 reads at a position in a bam file of 1 M reads is not the same thing as having 10 reads at the same position in a 100 M reads bam file. So I would advise to compute read depth by sample, normalize it by the library size of each sample and then to average over all samples if needed.
Then for the difference between merging bam and computing depth or computing depth for each bam and then averaging, it may come from the fact that in the first case, you're adding all the reads together, so what you get at the end is the sum and not the average (or the mean).
Otherwise, your method seems good for calculating read depths over a set of samples and of genes.
Good luck!