samtools depth print out all positions

I am trying to use samtools depth (v1.4) with the -a option and a bed file listing the human chromosomes chr1-chr22, chrX, chrY, and chrM to print out the coverage at every position:

cat GRCh38.karyo.bed | awk '{print $3}' | datamash sum 1 3088286401  I would like to know how to run samtools depth so that it produces 3,088,286,401 entries when run against a GRCh38 bam file: samtools depth -b$bedfile -a inputfile  I tried it for a few bam files that were aligned the same way, and I get differing number of entries: 3087003274 3087005666 3087007158 3087009435 3087009439 3087009621 3087009818 3087010065 3087010408 3087010477 3087010481 3087012115 3087013147 3087013186 3087013500 3087149616  Is there a special flag in samtools depth so that it reports all entries from the bed file? If samtools depth is not the best tool for this, what would be the equivalent with sambamba depth base? sambamba depth base --min-coverage=0 --regionsbedfile $inputfile  Any other options? • Out of interest, does -aa at least give the same number of bases each time? -a -a (or -aa) output absolutely all positions, including unused ref. sequences There is also an old but relevant discussion about -a and BED behaviour on github/samtools/#374 – Sam Nicholls Jun 7 '17 at 9:59 • I tried with -a and with -aa, and given a bed file, it seems to produce the same number of entries. – 719016 Jun 7 '17 at 12:35 • Two samtools commits from May 2016; Added 30ish new tests for depth/mpileup -a and -aa and Depth/Mpileup -a fixes. Fixes #374 are somewhat enlightening. My gut feeling is that although possible, using -a/-aa with a BED file is not entirely supported with samtools depth (i.e. has potentially undefined behaviour) and perhaps you would want to investigate another tool as you've suggested? – Sam Nicholls Jun 7 '17 at 13:43 • I dont work with human genomes, so excuse my naive question...why do you need to provide a bed file if you want every position? – Nathan S. Watson-Haigh Jun 7 '17 at 14:07 • If every reference defined in the BAM file is touched by at least one read, then -a and -aa will produce the same output. The second option, -aa is useful if you want to produce zero-depth information for sequences that are completely uncovered. – gringer Jun 7 '17 at 22:24 2 Answers You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file. bedtools genomecov -d -ibam$inputfile > "${inputfile}.genomecov"  You can also provide a BED file if you just want to calculate in the target region. Just replace the -a option with -aa: samtools depth -b$bedfile -aa \$inputfile


I see that you're using the GRCh38 human reference genome build, which includes alternate scaffolds that represent a wider variety of genomic variation in the human genome. This was quite a substantial change from previous reference genome builds, and many bioinformatics tools still haven't been updated to accept the GRCh38 reference.

The discrepancy in entry size is probably because there are alternate contigs that are uncovered by the reads that have been mapped. This is expected, but can cause a bit of confusion.

• Yes, the reason why I use GRCh38 bed file is to have the depth values for only the karyotypic chromosomes, rather than the scaffolds and decoys. The -aa option gave the same output as the -a option. – 719016 Jun 8 '17 at 8:45