samtools depth reports all positions though I don't indicate -a or -aa
0
0
Entering edit mode
4 weeks ago
vaushev ▴ 10

I am using a small test bam file; samtools flagstat does not show any issues:

$> samtools flagstat mytestbam.bam 20740 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 20740 + 0 mapped (100.00% : N/A) 20740 + 0 paired in sequencing 10260 + 0 read1 10480 + 0 read2 20118 + 0 properly paired (97.00% : N/A) 20303 + 0 with itself and mate mapped 437 + 0 singletons (2.11% : N/A) 185 + 0 with mate mapped to a different chr 185 + 0 with mate mapped to a different chr (mapQ>=5)  Now, I am trying to run samtools depth - and it reports all positions including those without any reads: $> samtools depth -r chr17:41196231-41196240 mytestbam.bam
chr17   41196231        0
chr17   41196232        0
chr17   41196233        0
chr17   41196234        0
chr17   41196235        0
chr17   41196236        0
chr17   41196237        0
chr17   41196238        1
chr17   41196239        1
chr17   41196240        1


I thought it should be the case only if I use -a or -aa flag. What am I doing wrong here, and how to make it report only positions with actual reads?

P.S. It's samtools 1.12, Using htslib 1.12

samtools • 177 views
0
Entering edit mode

I think there is no way to avoid position without reads, but I may be wrong. However, you can add an awk command to show only those positions with any reads samtools depth your_file.bam | awk '\$3 != 0'

0
Entering edit mode

Thank you Lila, this is a good workaround - but I still would like to understand why samtools behaves not the way it is described in the man. It clearly says that output all positions (including zero depth) is enabled by the -a flag.

0
Entering edit mode

I don't know why, but the manual said " -a option **may** sometimes operate... so as it said may... the answer is there, not always. But as I said, I can't tell. You always can try to contact the developers, in my experience, they are always very supportive. Good luck!