<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>

CC BY-SA 4.0
2020-08-09T02:33:37.840Z
2020-08-09T02:33:37.840Z
Mehrdad Bakhtiari
<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>

<sam><pysam>
Find all the bases for given reference position
# 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.


---

