Access base aligned to particular reference position

Question

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.

Answer 1

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.

Answer 2

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

Answer 3

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