BamOnBed

bob.py

import pysam
import os

t = "./865_rep1_treat.bam"
tb = "865_rep1_peaks.bed"

sorted_prefix = lambda path: path + ".sorted"
sorted_path = lambda path: sorted_prefix(path) + ".bam"


def make_index(in_bam, max_mem = 500000000):
    bai_path = lambda path: sorted_path(path) + ".bai"
    
    if not os.path.exists(sorted_path(in_bam)):
        pysam.sort("-m",str(max_mem),in_bam, sorted_prefix(in_bam))
    print "sorted"
        
    if not os.path.exists(bai_path(in_bam)):
        pysam.index(sorted_path(in_bam))
    print "indexed"
    return bai_path(in_bam)
               
result = make_index(t)

samfile = pysam.Samfile(sorted_path(t), "rb" )
with open(tb) as tbf:
    for i in tbf:
        chrom, peak_start, peak_end, = i.split("\t")[0:3]
        print [(i.pos,i.n) for i in samfile.pileup(chrom, int(peak_start), int(peak_end))]
samfile.close()