<p>By using <code>-h</code> in the <code>samtools view</code> command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the <code>-h</code> word count, but not the mpileup count. Try piping it through an additional <code>samtools view</code> (i.e. without <code>-h</code>) and see if the counts change:</p>
<pre><code>...
samtools view -h -f32 -  | \
samtools view | wc -l
</code></pre>
<p>Also, <code>samtools mpileup</code> by default only considers high-quality bases and concordant reads. Try adding a <code>-A</code> to your mpileup line (which stops anomalous read pairs from being discarded):</p>
<pre><code>...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l
</code></pre>
<p>Whether or not this is actually a good idea will be dependent on what you want to get out of the analysis, and what the downstream programs / analyses are expecting.</p>

CC BY-SA 3.0
2017-05-25T07:45:36.727Z
Bingo!  Adding `-A, --count-orphans     do not discard anomalous read pairs` did the trick!
719016
CC BY-SA 4.0
2024-02-29T15:05:52.657Z
i do have a follow up question. Is using the -A option a good idea? I have same problem as OP (i.e. empty mpileup), but can't find any useful info. of course with -A it seems to work but...how reliable will all downstream analysis be?
gbt
CC BY-SA 4.0
2024-02-29T20:37:04.653Z
As with most questions of this nature, there's not a universal answer. If you want a useful answer, please post this as a separate question with enough context to explain what the eventual end-goal is, and what that downstream analysis is.
gringer
CC BY-SA 4.0
2017-05-24T22:35:24.987Z
2024-02-29T20:34:37.437Z
2024-02-29T20:34:37.437Z
gringer
<p>I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.</p>

<p>I am heavily filtering reads that are not paired and of a certain orientation.
Applying the filtering, I get a few thousands of reads:</p>

<pre><code>samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l
</code></pre>

<p>Gives me <code>89502</code> reads.</p>

<p>If I then pipe this into <code>samtools mpileup</code>, I get no results:</p>

<pre><code>samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l
</code></pre>

<p>Returns 0.</p>

<p>I tried different combinations of filtering, and when I do both <code>-f 16</code> and <code>-f 32</code> returns empty, but if I do either of those, then it works:</p>

<pre><code>samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l
</code></pre>

<p>Returns <code>1056</code>.</p>

<p>Any ideas why? My thinking was that it would work with <code>--excl-flags 0</code>.</p>

<p>EDIT: substituting <code>mpileup</code> for <code>depth</code> does work, and prints out each position and the depth as expected.</p>

<p>EDIT2: adding <code>-q 0</code> to <code>mpileup</code> gives the same empty result.</p>

<p>Thanks in advance</p>

CC BY-SA 3.0
2017-05-24T16:52:07.737Z
Hi, could you add some details of how you got the synthetic reads and how you have aligned the to the reference?
Kamil S Jaron
CC BY-SA 3.0
2017-05-24T17:26:51.817Z
Are you aware that `-f 16` combined with `-f 32` (btw, you can just `-f 48` and skip the piping) are filtering for discordant reads? What happens if you add `-q 0` to mpileup?
Devon Ryan
CC BY-SA 3.0
2017-05-24T23:02:44.933Z
Since `depth` and `mpileup` share a fair bit of code, I have to wonder if you simply have no variants there. What happens if you specify `-a` and filter out the positions with no coverage?
Devon Ryan
CC BY-SA 3.0
2017-05-25T19:01:56.817Z
BTW, I suspect this will be a case where you'll need to post the BAM file somewhere. There are a few (currently commented out) `printf()`s in the samtools code that could probably be used to track down what's actually happening here.
Devon Ryan
CC BY-SA 3.0
2017-05-26T10:29:12.747Z
Devon already mentioned it in passing but it’s too important to get lost: Samtools flags are bit masks that you need to combine! Don’t call samtools repeatedly; call it once, with the combined bit flag. You can use the [BROAD online tool](https://broadinstitute.github.io/picard/explain-flags.html) to compute the combined flag if you don’t want to do it manually.
Konrad Rudolph
<sam><samtools><mpileup>
samtools mpileup empty when filtering out flags
# samtools mpileup empty when filtering out flags

(Asked by: username_0 on 2017-05-24)

I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.
I am heavily filtering reads that are not paired and of a certain orientation.
Applying the filtering, I get a few thousands of reads:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l

Gives me 89502 reads.
If I then pipe this into samtools mpileup, I get no results:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 0.
I tried different combinations of filtering, and when I do both -f 16 and -f 32 returns empty, but if I do either of those, then it works:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 1056.
Any ideas why? My thinking was that it would work with --excl-flags 0.
EDIT: substituting mpileup for depth does work, and prints out each position and the depth as expected.
EDIT2: adding -q 0 to mpileup gives the same empty result.
Thanks in advance

**Tags:** sam, samtools, mpileup
> **Comments:**
> **username_1**: Hi, could you add some details of how you got the synthetic reads and how you have aligned the to the reference?  
> **username_2**: Are you aware that `-f 16` combined with `-f 32` (btw, you can just `-f 48` and skip the piping) are filtering for discordant reads? What happens if you add `-q 0` to mpileup?  
> **username_2**: Since `depth` and `mpileup` share a fair bit of code, I have to wonder if you simply have no variants there. What happens if you specify `-a` and filter out the positions with no coverage?  
> **username_2**: BTW, I suspect this will be a case where you'll need to post the BAM file somewhere. There are a few (currently commented out) `printf()`s in the samtools code that could probably be used to track down what's actually happening here.  
> **username_3**: Devon already mentioned it in passing but it’s too important to get lost: Samtools flags are bit masks that you need to combine! Don’t call samtools repeatedly; call it once, with the combined bit flag. You can use the [BROAD online tool](https://broadinstitute.github.io/picard/explain-flags.html) to compute the combined flag if you don’t want to do it manually.  

---

## 1 answers

### Answer by username_4 on 2017-05-24. Score: 2 | Accepted answer

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts change:
...
samtools view -h -f32 -  | \
samtools view | wc -l

Also, samtools mpileup by default only considers high-quality bases and concordant reads. Try adding a -A to your mpileup line (which stops anomalous read pairs from being discarded):
...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l

Whether or not this is actually a good idea will be dependent on what you want to get out of the analysis, and what the downstream programs / analyses are expecting.

> **Comments:**
> **username_0**: Bingo!  Adding `-A, --count-orphans     do not discard anomalous read pairs` did the trick!  
> **username_5**: i do have a follow up question. Is using the -A option a good idea? I have same problem as OP (i.e. empty mpileup), but can't find any useful info. of course with -A it seems to work but...how reliable will all downstream analysis be?  
> **username_4**: As with most questions of this nature, there's not a universal answer. If you want a useful answer, please post this as a separate question with enough context to explain what the eventual end-goal is, and what that downstream analysis is.  

---

