Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-25T07:45:36.727Z
Answers.Text: Bingo!  Adding `-A, --count-orphans     do not discard anomalous read pairs` did the trick!
Answers.UserDisplayName: 719016
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-29T15:05:52.657Z
Answers.Text: 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?
Answers.UserDisplayName: gbt
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-02-29T20:37:04.653Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-05-24T22:35:24.987Z
Answers.LastActivityDate: 2024-02-29T20:34:37.437Z
Answers.LastEditDate: 2024-02-29T20:34:37.437Z
Answers.UserDisplayName: gringer
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-24T16:52:07.737Z
Comments.Text: Hi, could you add some details of how you got the synthetic reads and how you have aligned the to the reference?
Comments.UserDisplayName: Kamil S Jaron
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-24T17:26:51.817Z
Comments.Text: 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?
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-24T23:02:44.933Z
Comments.Text: 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?
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-25T19:01:56.817Z
Comments.Text: 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.
Comments.UserDisplayName: Devon Ryan
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-26T10:29:12.747Z
Comments.Text: 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.
Comments.UserDisplayName: Konrad Rudolph
Tags: <sam><samtools><mpileup>
Title: samtools mpileup empty when filtering out flags
ThreadText: # 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**: username_2 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.  
  
  ---
  
  