#!/usr/bin/env python2

import sys,subprocess
from guppy import hpy
from subprocess import Popen

## Reads a bam file from stdin and compares it to a modified bam file. for each read, the location of original and modified BAM are compared and only reads mapping to the same position are written to stdout
##samtools calmd $originalbam $refseq | samtools view -L $bedfile -h -q 30 - | filter_sam_startpos_dict.py $modifiedbam | samtools view -bh - > $modified.filter.bam


bamfile=sys.argv[1]

start_dict={}

cmd="samtools view -q 20 %s" %(bamfile)
	
process = Popen(cmd, shell=True, stdout=subprocess.PIPE,stderr=subprocess.PIPE)

# first goes through the modified bam (argument 1) and makes a dictionary of read names, chroms and positions
for l in iter(process.stdout.readline,''):
	split=l.split()
	md=l.split('MD:Z:')[1].split('\t')[0]
	if md.find('^')<0:
		start_dict[split[0].replace("mod_",'')]=(split[2],split[3])
	



# next goes through the original bam (stdin) and checks if the chrom and star match the dictionary. Throws out read if md tag has insertion. mismatches variable seems to be unused.
for line in sys.stdin:
	if line[0]=='@':
		print line.rstrip()
		continue
	split=line.split()
	read_ID=split[0]
	c=split[2]
	bp=split[3]

	found=0

	if start_dict.has_key(read_ID):
		t=start_dict[read_ID]
		if c==t[0] and bp ==t[1]:
			if line.find('NM:i:')>=0 and line.find('MD:Z:'):
				mismatches=int(line.split('NM:i:')[1].split('\t')[0])
	                	md=line.split('MD:Z:')[1].split('\t')[0]
	        	        if md.find('^')<0:
        	        		found=1

	if found:
		print line.rstrip()

#h = hpy()
#f=open('memory.info','w')
#f.write(str(h.heap()))
#f.close()
