Using pybedtools to report a BAM read's allele when it overlaps a SNP

get_bam_allele.py

#!/usr/bin/env python

"""
get_bam_allele.py
Usage: python get_bam_allele.py [BAM] [BED]
"""

from pybedtools import BedTool
import sys

# Create "BedTools" for the BAM and BED file
# that are passed in from the command line.
reads = BedTool(sys.argv[1])
snps  = BedTool(sys.argv[2])

# Loop through each read in the BAM file and
# report the base for each SNP that it overlaps
for read in reads:
    for snp in snps.all_hits(read):
        # print the read name (idx=0), snp name and the allele
        # in the read at the SNP site (idx=9 is the BAM sequence)
        print read[0], snp.name, read[9][snp.start - read.start]