Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T08:45:00.740Z
Answers.Text: 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.
Answers.UserDisplayName: 719016
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-07T22:34:32.847Z
Answers.LastActivityDate: 2017-06-07T22:34:32.847Z
Answers.UserDisplayName: gringer
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-07T16:20:19.863Z
Answers.LastActivityDate: 2017-06-07T16:20:19.863Z
Answers.UserDisplayName: morgantaschuk
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T09:59:26.187Z
Comments.Text: 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)
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T12:35:22.353Z
Comments.Text: I tried with `-a` and with `-aa`, and given a bed file, it seems to produce the same number of entries.
Comments.UserDisplayName: 719016
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T13:43:32.487Z
Comments.Text: 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?
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T14:07:18.563Z
Comments.Text: 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?
Comments.UserDisplayName: Nathan S. Watson-Haigh
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-07T22:24:57.373Z
Comments.Text: 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.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T08:45:31.897Z
Comments.Text: @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.
Comments.UserDisplayName: 719016
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-08T21:07:18.540Z
Comments.Text: ... except it appears that that there's a bug in samtools around this (i.e. I just noticed Sam Nicholls' comment).
Comments.UserDisplayName: gringer
Tags: <sam><samtools><sambamba>
Title: samtools depth print out all positions
ThreadText: # 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**: @username_3 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 username_1' comment).  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_4 on 2017-06-07. Score: 6 | Accepted answer
  
  You might try using [bedtools genomecov](http://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html) 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](https://www.quora.com/What-are-the-differences-between-GRCh38-and-GRCh37), 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.  
  
  ---
  
  