Answers.Body: <p><code>pileupcolumn.n</code> shows the total number of the reads that cover this basepair (in this case 24793). However, <code>pileupcolumn.pileups</code> iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting <code>min_base_quality</code> parameter when you call pileup.</p>
  <p>To have a similar filter applied to both coverage and reads, you can use <code>piliecolumn.get_num_aligned()</code> to get the number of segments and <code>pileupcolumn.pileups</code> to get the segments themselves.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-09T02:33:37.840Z
Answers.LastActivityDate: 2020-08-09T02:33:37.840Z
Answers.UserDisplayName: Mehrdad Bakhtiari
Body: <p>Hi im looking for the aligned bases in the reads for a given reference position. im using the following script from the pysam documentataion.
  I adjusted it to find the specified position. in this case 24793.</p>
  <pre><code>import pysam
  
  samfile = pysam.AlignmentFile(&quot;generated_alignment_sorted.bam&quot;, &quot;rb&quot; )
  for pileupcolumn in samfile.pileup(&quot;chr05_modified.copy0&quot;, 10, 52000000):
      if pileupcolumn.pos == 24793:
          print (&quot;\ncoverage at base %s = %s&quot; %
                 (pileupcolumn.pos, pileupcolumn.n))
          for pileupread in pileupcolumn.pileups:
                  print ('\tbase in read %s = %s' %
                        (pileupread.alignment.query_name.split(';')[0],
                         pileupread.alignment.query_sequence[pileupread.query_position]))
  
  samfile.close()
  </code></pre>
  <p>This outputs:</p>
  <pre><code>coverage at base 24793 = 6
          base in read m419941/6207/CCS Read=419941 = C
  </code></pre>
  <p>For all the position I tried it doesnt print out the same number of bases as it says that the coverage is.
  shouldnt this print out 6 times: base in read m419941/6207/CCS Read=419941 = C for different reads?</p>
  <p>Why is this?</p>
  
Tags: <sam><pysam>
Title: Find all the bases for given reference position
ThreadText: # Find all the bases for given reference position
  
  (Asked by: username_0 on 2020-07-31)
  
  Hi im looking for the aligned bases in the reads for a given reference position. im using the following script from the pysam documentataion.
  I adjusted it to find the specified position. in this case 24793.
  
  ```
  import pysam
  
  samfile = pysam.AlignmentFile("generated_alignment_sorted.bam", "rb" )
  for pileupcolumn in samfile.pileup("chr05_modified.copy0", 10, 52000000):
      if pileupcolumn.pos == 24793:
          print ("\ncoverage at base %s = %s" %
                 (pileupcolumn.pos, pileupcolumn.n))
          for pileupread in pileupcolumn.pileups:
                  print ('\tbase in read %s = %s' %
                        (pileupread.alignment.query_name.split(';')[0],
                         pileupread.alignment.query_sequence[pileupread.query_position]))
  
  samfile.close()
  
  ```
  
  This outputs:
  
  ```
  coverage at base 24793 = 6
          base in read m419941/6207/CCS Read=419941 = C
  
  ```
  
  For all the position I tried it doesnt print out the same number of bases as it says that the coverage is.
  shouldnt this print out 6 times: base in read m419941/6207/CCS Read=419941 = C for different reads?
  
  Why is this?
  
  **Tags:** sam, pysam
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-08-09. Score: 1 | Accepted answer
  
  `pileupcolumn.n` shows the total number of the reads that cover this basepair (in this case 24793). However, `pileupcolumn.pileups` iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting `min_base_quality` parameter when you call pileup.
  
  To have a similar filter applied to both coverage and reads, you can use `piliecolumn.get_num_aligned()` to get the number of segments and `pileupcolumn.pileups` to get the segments themselves.
  
  
  ---
  
  