<p>You can use the <a href="https://pysam.readthedocs.io/en/latest/api.html#pysam.AlignedSegment.get_reference_positions" rel="nofollow noreferrer"><code>get_reference_positions()</code></a> function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.</p>

<pre><code>for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...
</code></pre>

<p>That's not exactly correct syntax, but you get the idea.</p>

CC BY-SA 4.0
2019-04-08T14:57:18.830Z
Thanks, this is helpful. When there are no insertions (`I` blocks) in the CIGAR, I can use `idx` and `idx2` to index directly in to `record.query_sequence` and get the base at the correct position. However, reads with insertions need additional operations to line up the `get_reference_positions()` indices with the `query_sequence` indices. So unless I'm missing something, it looks like there's a way to do this without computing a pileup...just not a simple way. :-)
Daniel Standage
CC BY-SA 4.0
2019-04-08T14:58:49.607Z
I think you summarized it nicely, possible but not simple.
Devon Ryan
CC BY-SA 4.0
2019-04-07T07:26:36.660Z
2019-04-08T14:21:13.610Z
2019-04-08T14:21:13.610Z
Devon Ryan
<p>I had the same problem as you, and so I wrote my own code and sharing it here. No guarantee that it's bug-free ;-)</p>
<p>This code takes advantage of the cigar string. It returns None if the base has been deleted or clipped, or if the bam file did not store the sequence.</p>
<pre><code>from typing import Optional
import pysam



def find_base_in_alignment(alignment: pysam.AlignedSegment,
                           pos: int, 
                           bam_stores_revcomp: bool = False) \ 
        -&gt; Optional[str]:
    idx_q = 0
    idx_r = pos - alignment.reference_start
    if bam_stores_revcomp:
        seq = alignment.query_sequence
    else:
        seq = alignment.get_forward_sequence()
    
    if seq is None:
        return None
    
    for op, l in alignment.cigartuples:
        ref_consumed = op in {0, 2, 3, 7, 8}
        query_consumed = op in {0, 1, 4, 7, 8}
        
        if ref_consumed:
            idx_r -= l
        if query_consumed:
            idx_q += l
        
        if idx_r &lt; 0:
            if query_consumed:
                # base is in query between idx_q-l , idx_q
                base = seq[idx_q + idx_r - 1]
                return base
            else:
                # position has been deleted
                return None
</code></pre>
<p>Let's say you want to genotype all reads by looking at the base at chromosome <code>chrom</code> and position <code>pos</code>:</p>
<pre><code>bamfile = pysam.AlignmentFile(filename, &quot;r&quot; )
for alignment in bamfile.fetch(chrom, pos-1, pos+1):
    base = find_base_in_alignment(alignment, pos)
    # ... do something with it
</code></pre>
<p>Note that the parameter <code>bam_stores_revcomp</code> is important. Some aligners will store the original query sequence in the bam, and other aligners (like minimap2) will store the forward sequence (i.e. the reverse complement whenever a read was aligned to the reverse strand). If you use the former, set <code>bam_stores_revcomp</code> to <code>False</code>, otherwise set it to <code>True</code></p>

CC BY-SA 4.0
2021-06-26T01:40:04.967Z
Hi, thank you for your code! I think it should be `if idx_r <= 0:` , then your code would report the exact same data as IGV, excluding insertions. Thank you!
Roberto
CC BY-SA 4.0
2021-08-29T01:17:50.040Z
this should be the accepted answer.
Maximilian Press
CC BY-SA 4.0
2021-05-26T16:10:40.643Z
2021-05-27T09:13:49.007Z
2021-05-27T09:13:49.007Z
RHS
<p>This is not [yet] a complete answer, but hopefully it helps get you on the right path.</p>

<p>I wrote <a href="https://gitlab.com/gringer/bioinfscripts/blob/master/sam2mapped.pl" rel="nofollow noreferrer">some code</a> to retrieve the subsequence from reads in a BAM file that are actually mapped to the reference, which required calculating mapped base locations. Here's a subsection of the code that does this (slightly paraphrased):</p>

<pre><code>  my <span class="math-container">$refPos = $</span>F[3];
  my <span class="math-container">$cigar = $</span>F[5];
  <span class="math-container">$cigar =~ s/[0-9]S$</span>//;
  my <span class="math-container">$seq = $</span>F[9];
  my <span class="math-container">$qual = $</span>F[10];
  my <span class="math-container">$startTrim = 0;
  my $</span>matchLen = 0;
  while(<span class="math-container">$cigar =~ s/^([0-9]+)([MIDNSHP=X])//){
    my $subLen = $1;
    my $op = $2;
    if($op eq "S"){
        $seq = substr($seq, $subLen);
        $qual = substr($qual, $subLen);
      }
    }
    if($</span>op =~ /[M=XI]/){
      <span class="math-container">$matchLen += $</span>subLen;
    }
  }
  <span class="math-container">$seq = substr($</span>seq, 0, <span class="math-container">$matchLen + $</span>addSeq);
  <span class="math-container">$qual = substr($</span>qual, 0, <span class="math-container">$matchLen + $</span>addSeq);
</code></pre>

<p>So, some things to take note of:</p>

<ul>
<li>The entire CIGAR string (excluding terminal soft clips) should be parsed up to the required base location</li>
<li>Any soft clips at the start of the sequence need to be accounted for</li>
<li>Only update the read location counter if there's a match, mismatch, or insertion</li>
</ul>

CC BY-SA 4.0
2019-04-09T09:58:48.090Z
2019-04-09T09:58:48.090Z
gringer
<p><strong>The short version</strong>: If I have a SAM record, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?</p>

<p><strong>The long version</strong>: I'm using pysam to write some genotyping code. I have a BAM file with reads aligned to some amplicons. I'd now like to genotype a couple of positions of interest in each read while maintaining the phase information. For example, at amplicon 1, if read 1 has <code>A</code> for SNP1 and <code>T</code> at SNP2, then I want to report an <code>A,T</code> genotype for read 1 at amplicon 1. </p>

<p>The typical way to genotype SNPs of interest is to compute a pileup and look at all of the aligned bases position-by-position. The pysam API includes the following example of how to work with pileups.</p>

<pre class="lang-py prettyprint-override"><code>import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()
</code></pre>

<p>However, I lose the phasing information if I analyze the data position-by-position. Is there any way to flip this example so that I'm iterating over reads in the outer loop, and then over positions in the inner loop? In other words, <strong>if I'm iterating over reads, is there any simple way to retrieve the base aligned to a particular reference position <em>without</em> computing a pileup</strong>?</p>

<p>As a potential solution, I've though about storing genotype information in a nested dictionary indexed first by read ID and then by position, but depending on coverage that could be incredibly memory intensive.</p>

CC BY-SA 4.0
2023-09-30T12:39:31.463Z
You can use `get_aligned_pairs()`; you'll get a list of tuples of `(position_in_query_sequence, position_in_reference_sequence)`
Fravadona
CC BY-SA 4.0
2023-10-02T17:50:06.830Z
That's promising! It's been a hot second since I've worked on this code, but it looks like it could indeed allow me to keep phase information on a read-by-read basis. I'll try it out when I can.
Daniel Standage
<sam><pysam>
Access base aligned to particular reference position
# Access base aligned to particular reference position

(Asked by: username_0 on 2019-04-05)

The short version: If I have a SAM record, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?
The long version: I'm using pysam to write some genotyping code. I have a BAM file with reads aligned to some amplicons. I'd now like to genotype a couple of positions of interest in each read while maintaining the phase information. For example, at amplicon 1, if read 1 has A for SNP1 and T at SNP2, then I want to report an A,T genotype for read 1 at amplicon 1. 
The typical way to genotype SNPs of interest is to compute a pileup and look at all of the aligned bases position-by-position. The pysam API includes the following example of how to work with pileups.
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

However, I lose the phasing information if I analyze the data position-by-position. Is there any way to flip this example so that I'm iterating over reads in the outer loop, and then over positions in the inner loop? In other words, if I'm iterating over reads, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?
As a potential solution, I've though about storing genotype information in a nested dictionary indexed first by read ID and then by position, but depending on coverage that could be incredibly memory intensive.

**Tags:** sam, pysam
> **Comments:**
> **username_1**: You can use `get_aligned_pairs()`; you'll get a list of tuples of `(position_in_query_sequence, position_in_reference_sequence)`  
> **username_0**: That's promising! It's been a hot second since I've worked on this code, but it looks like it could indeed allow me to keep phase information on a read-by-read basis. I'll try it out when I can.  

---

## 3 answers

### Answer by username_2 on 2019-04-07. Score: 4

You can use the get_reference_positions() function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.
for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...

That's not exactly correct syntax, but you get the idea.

> **Comments:**
> **username_0**: Thanks, this is helpful. When there are no insertions (`I` blocks) in the CIGAR, I can use `idx` and `idx2` to index directly in to `record.query_sequence` and get the base at the correct position. However, reads with insertions need additional operations to line up the `get_reference_positions()` indices with the `query_sequence` indices. So unless I'm missing something, it looks like there's a way to do this without computing a pileup...just not a simple way. :-)  
> **username_2**: I think you summarized it nicely, possible but not simple.  

---

### Answer by username_3 on 2021-05-26. Score: 3

I had the same problem as you, and so I wrote my own code and sharing it here. No guarantee that it's bug-free ;-)
This code takes advantage of the cigar string. It returns None if the base has been deleted or clipped, or if the bam file did not store the sequence.
from typing import Optional
import pysam



def find_base_in_alignment(alignment: pysam.AlignedSegment,
                           pos: int, 
                           bam_stores_revcomp: bool = False) \ 
        -> Optional[str]:
    idx_q = 0
    idx_r = pos - alignment.reference_start
    if bam_stores_revcomp:
        seq = alignment.query_sequence
    else:
        seq = alignment.get_forward_sequence()
    
    if seq is None:
        return None
    
    for op, l in alignment.cigartuples:
        ref_consumed = op in {0, 2, 3, 7, 8}
        query_consumed = op in {0, 1, 4, 7, 8}
        
        if ref_consumed:
            idx_r -= l
        if query_consumed:
            idx_q += l
        
        if idx_r < 0:
            if query_consumed:
                # base is in query between idx_q-l , idx_q
                base = seq[idx_q + idx_r - 1]
                return base
            else:
                # position has been deleted
                return None

Let's say you want to genotype all reads by looking at the base at chromosome chrom and position pos:
bamfile = pysam.AlignmentFile(filename, "r" )
for alignment in bamfile.fetch(chrom, pos-1, pos+1):
    base = find_base_in_alignment(alignment, pos)
    # ... do something with it

Note that the parameter bam_stores_revcomp is important. Some aligners will store the original query sequence in the bam, and other aligners (like minimap2) will store the forward sequence (i.e. the reverse complement whenever a read was aligned to the reverse strand). If you use the former, set bam_stores_revcomp to False, otherwise set it to True

> **Comments:**
> **username_4**: Hi, thank you for your code! I think it should be `if idx_r <= 0:` , then your code would report the exact same data as IGV, excluding insertions. Thank you!  
> **username_5**: this should be the accepted answer.  

---

### Answer by username_6 on 2019-04-09. Score: 1

This is not [yet] a complete answer, but hopefully it helps get you on the right path.
I wrote some code to retrieve the subsequence from reads in a BAM file that are actually mapped to the reference, which required calculating mapped base locations. Here's a subsection of the code that does this (slightly paraphrased):
  my $refPos = $F[3];
  my $cigar = $F[5];
  $cigar =~ s/[0-9]S$//;
  my $seq = $F[9];
  my $qual = $F[10];
  my $startTrim = 0;
  my $matchLen = 0;
  while($cigar =~ s/^([0-9]+)([MIDNSHP=X])//){
    my $subLen = $1;
    my $op = $2;
    if($op eq "S"){
        $seq = substr($seq, $subLen);
        $qual = substr($qual, $subLen);
      }
    }
    if($op =~ /[M=XI]/){
      $matchLen += $subLen;
    }
  }
  $seq = substr($seq, 0, $matchLen + $addSeq);
  $qual = substr($qual, 0, $matchLen + $addSeq);

So, some things to take note of:

The entire CIGAR string (excluding terminal soft clips) should be parsed up to the required base location
Any soft clips at the start of the sequence need to be accounted for
Only update the read location counter if there's a match, mismatch, or insertion


---

