text: user96: 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.
  
  user73: 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
  
  
  user12344: 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
  
  user77: 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.
  
original_text: 96: <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>
   73: <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>
   12344: <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>
   77: <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>
  