Answers.Body: <p>The two entries in the sam file represent mate pairs for the given ID.  You can tell the difference based on the sam flags.  <a href="http://broadinstitute.github.io/picard/explain-flags.html" rel="nofollow noreferrer">Picard</a> has a nice tool to determine the meaning of the flags.</p>
  
  <p>You are missing coverage within the first read.  I think that it is related to the low quality bases.  samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-06T00:53:54.123Z
Answers.Text: Thanks for the tip about distinguishing a read and its mate by QNAME _and_ flag—that's helpful. The read with flag `83`, though, has `baseQ > 13` at 150 positions, yet appears in the pileup columns for only 121 positions.
Answers.UserDisplayName: Randoms
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-06T02:17:15.543Z
Answers.Text: Ah but **BAQ** ("Base Alignment Quality"), not just regular base quality, also has to be above 13. Comparing `samtools mpileup` and `samtools mpileup -Q0` shows that that's why reads are being filtered out.
Answers.UserDisplayName: Randoms
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-05T22:30:21.297Z
Answers.LastActivityDate: 2018-07-05T22:30:21.297Z
Answers.UserDisplayName: Bioathlete
Body: <p><a href="https://www.biostars.org/p/324857/" rel="nofollow noreferrer">Cross-posted</a> on Biostars (with no answers currently). I hope that's OK. </p>
  
  <hr>
  
  <p>I thought that <strong><em>iff</em></strong> a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case. </p>
  
  <p>This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.</p>
  
  <p>What's wrong in my understanding?</p>
  
  <hr>
  
  <p>Updated code (and updated output) moved to <a href="https://gist.github.com/jaLasser/15ce9d5f6b91854fe672b5a830ed4f5f" rel="nofollow noreferrer">Gist</a>.</p>
  
  <hr>
  
  <p><strong>Edit</strong>: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.</p>
  
  <pre><code>E00489:44:HNNVYCCXX:1:2120:26524:19135  163     5       345248  0       48M102S =       345446  350     TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG  AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)&lt;)&lt;--)&lt;)--&lt;-FAJFJFFFF-&lt;  XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M   MD:Z:48 RG:Z:HNNVYCCXX.1       NM:i:0   MQ:i:60 AS:i:48 XS:i:48
  E00489:44:HNNVYCCXX:1:2120:26524:19135  83      5       345446  60      6M2D144M        =       345248  -350    TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG  JJAF7-FJJJFJFF-AFJF7--&lt;A--&lt;F&lt;7-7--FA--F-J--A-7--7-FJFJ7JFJJF&lt;JA77-7----7----&lt;F7-&lt;JJF-&lt;&lt;FAJJJJA-7-&lt;JAJJ&lt;7FFJJ&lt;JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MC:Z:48M102S    MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1        NM:i:7  MQ:i:0  AS:i:119        XS:i:69
  </code></pre>
  
  <p>I guess this explains why the read appears in more than <code>150</code> pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in <code>198</code> = (<code>48</code> + <code>6</code> + <code>144</code>) pileup columns, at each position where one of them maps?</p>
  
  <hr>
  
  <p><strong>Update:</strong> I re-ran this script, using @Bioathlete's suggestion of using the flags to distinguish the read and its mate. In case <code>pileup()</code> was filtering reads with low base quality, I also calculated how many positions the read has <code>baseQ &gt;= 13</code>. It's <code>150</code> positions, but this read still appears in pileup columns for only <code>121</code> positions. </p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-07-05T19:06:45.167Z
Comments.Text: Can you include the raw sam lines for the read so we can determine if the issue is in pysam or your code.  That would give us truth.
Comments.UserDisplayName: Bioathlete
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-07-05T19:16:49.967Z
Comments.Text: @Bioathlete please see above
Comments.UserDisplayName: Randoms
Tags: <python><sam><phylogenetics><pysam>
Title: pysam pileup: what reads appear in the pileup?
ThreadText: # pysam pileup: what reads appear in the pileup?
  
  (Asked by: username_0 on 2018-07-05)
  
  [Cross-posted](https://www.biostars.org/p/324857/) on Biostars (with no answers currently). I hope that's OK.
  
  ---
  
  I thought that ***iff*** a read matches the reference at a position, it would appear in the pileup column for that position. But tracking a single read in pysam shows that that's not the case.
  
  This read appears in some pileup columns for positions to which it does not map, and does not appear in some pileup columns to which it does map.
  
  What's wrong in my understanding?
  
  ---
  
  Updated code (and updated output) moved to [Gist](https://gist.github.com/jaLasser/15ce9d5f6b91854fe672b5a830ed4f5f).
  
  ---
  
  **Edit**: Here's the raw SAM read. I wasn't separating the read from it's mate when I calculated the figures above.
  
  ```
  E00489:44:HNNVYCCXX:1:2120:26524:19135  163     5       345248  0       48M102S =       345446  350     TGTGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTGTGGGCGGGAGTGCGTGGGGGAGGGTGAGAGTGTGGATGTCGGTGGTCGCGGCTGCGTGAGAGAGGGGGTGTGTGGGGGGGGGGGGGGGGGGGGGTGTGGGTGGGCG  AAAFFJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJ-J-7-777-7-7-7--7-77-77F---7-7AFJ7----7-7FAA7-7---7----77-777------7-7)7-AAA-A)7FA)<)<--)<)--<-FAJFJFFFF-<  XA:Z:5,+345312,48M102S,0;5,+345554,48M102S,1;5,+345452,4S44M102S,0;5,+345192,48M102S,1; MC:Z:6M2D144M   MD:Z:48 RG:Z:HNNVYCCXX.1       NM:i:0   MQ:i:60 AS:i:48 XS:i:48
  E00489:44:HNNVYCCXX:1:2120:26524:19135  83      5       345446  60      6M2D144M        =       345248  -350    TGGGGATGTGTGTGTGTGTGTCGGATGATGTCCCTGGTTGTGTGTGGGGATGTGTGTGTGTGTGGGTTGATGGTCCTGGCTGTGTGTGTGTGTGGGTGTGCGTGTGTGGGTGTGTGTGTGTGTGTGTCGGATGATGTCCCTGGCTGTGTG  JJAF7-FJJJFJFF-AFJF7--<A--<F<7-7--FA--F-J--A-7--7-FJFJ7JFJJF<JA77-7----7----<F7-<JJF-<<FAJJJJA-7-<JAJJ<7FFJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFFFAA  MC:Z:48M102S    MD:Z:6^TG31C25C2A5T0C76 RG:Z:HNNVYCCXX.1        NM:i:7  MQ:i:0  AS:i:119        XS:i:69
  
  ```
  
  I guess this explains why the read appears in more than `150` pileup columns—in some of the columns, it's really its mate. But given the CIGAR strings for these paired reads, I would expect one of these mates to appear in `198` = (`48` + `6` + `144`) pileup columns, at each position where one of them maps?
  
  ---
  
  **Update:** I re-ran this script, using @username_1's suggestion of using the flags to distinguish the read and its mate. In case `pileup()` was filtering reads with low base quality, I also calculated how many positions the read has `baseQ >= 13`. It's `150` positions, but this read still appears in pileup columns for only `121` positions.
  
  **Tags:** python, sam, phylogenetics, pysam
  
  > **Comments:**
  > **username_1**: Can you include the raw sam lines for the read so we can determine if the issue is in pysam or your code.  That would give us truth.  
  > **username_0**: @username_1 please see above  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2018-07-05. Score: 2 | Accepted answer
  
  The two entries in the sam file represent mate pairs for the given ID. You can tell the difference based on the sam flags. [Picard](http://broadinstitute.github.io/picard/explain-flags.html) has a nice tool to determine the meaning of the flags.
  
  You are missing coverage within the first read. I think that it is related to the low quality bases. samtools mpileup filters bases with a quality below 13 and pysam just wraps the samtools C functions.
  
  > **Comments:**
  > **username_0**: Thanks for the tip about distinguishing a read and its mate by QNAME _and_ flag—that's helpful. The read with flag `83`, though, has `baseQ > 13` at 150 positions, yet appears in the pileup columns for only 121 positions.  
  > **username_0**: Ah but **BAQ** ("Base Alignment Quality"), not just regular base quality, also has to be above 13. Comparing `samtools mpileup` and `samtools mpileup -Q0` shows that that's why reads are being filtered out.  
  
  ---
  
  