<p>Just replace the <code>-a</code> option with <code>-aa</code>:</p>

<pre><code>samtools depth -b $bedfile -aa $inputfile
</code></pre>

<p>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 <a href="https://www.quora.com/What-are-the-differences-between-GRCh38-and-GRCh37" rel="nofollow noreferrer">previous reference genome builds</a>, and many bioinformatics tools still haven't been updated to accept the GRCh38 reference.</p>

<p>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.</p>

CC BY-SA 3.0
2017-06-08T08:45:00.740Z
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
CC BY-SA 3.0
2017-06-07T22:34:32.847Z
2017-06-07T22:34:32.847Z
gringer
<p>You might try using <a href="http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html" rel="noreferrer">bedtools genomecov</a> instead. If you provide the -d option, it reports the coverage at every position in the BAM file.</p>

<pre><code>bedtools genomecov -d -ibam $inputfile &gt; "${inputfile}.genomecov"
</code></pre>

<p>You can also provide a BED file if you just want to calculate in the target  region.</p>

CC BY-SA 3.0
2017-06-07T16:20:19.863Z
2017-06-07T16:20:19.863Z
morgantaschuk
<p>I am trying to use <code>samtools depth</code> (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:</p>

<pre><code>cat GRCh38.karyo.bed | awk '{print $3}' | datamash sum 1
3088286401
</code></pre>

<p>I would like to know how to run <code>samtools depth</code> so that it produces 3,088,286,401 entries when run against a GRCh38 bam file:</p>

<pre><code>samtools depth -b $bedfile -a $inputfile
</code></pre>

<p>I tried it for a few bam files that were aligned the same way, and I get differing number of entries:</p>

<pre><code>3087003274
3087005666
3087007158
3087009435
3087009439
3087009621
3087009818
3087010065
3087010408
3087010477
3087010481
3087012115
3087013147
3087013186
3087013500
3087149616
</code></pre>

<p>Is there a special flag in <code>samtools depth</code> so that it reports all entries from the bed file?</p>

<p>If <code>samtools depth</code> is not the best tool for this, what would be the equivalent with <code>sambamba depth base</code>?</p>

<pre><code>sambamba depth base --min-coverage=0 --regions $bedfile $inputfile
</code></pre>

<p>Any other options?</p>

CC BY-SA 3.0
2017-06-07T09:59:26.187Z
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](https://github.com/samtools/samtools/issues/374)
Sam Nicholls
CC BY-SA 3.0
2017-06-07T12:35:22.353Z
I tried with `-a` and with `-aa`, and given a bed file, it seems to produce the same number of entries.
719016
CC BY-SA 3.0
2017-06-07T13:43:32.487Z
Two `samtools` commits from May 2016; [Added 30ish new tests for depth/mpileup -a and -aa](https://github.com/jkbonfield/samtools/commit/b4a5902d2430a83ab63ce84be48a5acc52458693) and [Depth/Mpileup -a fixes. Fixes #374](https://github.com/samtools/samtools/commit/44290ddc3c8808cd82689d6f206c47fb112a59a1) 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
CC BY-SA 3.0
2017-06-07T14:07:18.563Z
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
CC BY-SA 3.0
2017-06-07T22:24:57.373Z
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
CC BY-SA 3.0
2017-06-08T08:45:31.897Z
@gringer 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.
719016
CC BY-SA 3.0
2017-06-08T21:07:18.540Z
... except it appears that that there's a bug in samtools around this (i.e. I just noticed Sam Nicholls' comment).
gringer
<sam><samtools><sambamba>
samtools depth print out all positions
# samtools depth print out all positions

(Asked by: username_0 on 2017-06-07)

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 --regions $bedfile $inputfile

Any other options?

**Tags:** sam, samtools, sambamba
> **Comments:**
> **username_1**: 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](https://github.com/samtools/samtools/issues/374)  
> **username_0**: I tried with `-a` and with `-aa`, and given a bed file, it seems to produce the same number of entries.  
> **username_1**: Two `samtools` commits from May 2016; [Added 30ish new tests for depth/mpileup -a and -aa](https://github.com/jkbonfield/samtools/commit/b4a5902d2430a83ab63ce84be48a5acc52458693) and [Depth/Mpileup -a fixes. Fixes #374](https://github.com/samtools/samtools/commit/44290ddc3c8808cd82689d6f206c47fb112a59a1) 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?  
> **username_2**: 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?  
> **username_3**: 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.  
> **username_0**: @gringer 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.  
> **username_3**: ... except it appears that that there's a bug in samtools around this (i.e. I just noticed Sam Nicholls' comment).  

---

## 2 answers

### Answer by username_4 on 2017-06-07. Score: 6 | Accepted answer

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.


---

### Answer by username_3 on 2017-06-07. Score: 2

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.

> **Comments:**
> **username_0**: 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.  

---

