#!/usr/bin/python

"""
   This script is an HTSeq alternatives.
   It is count reads in regions, but if there
   are overlapping regions the read is counted in both
   regions. Read also counted if partially overlap
"""

import sys
import gzip

class Exon:
	def __init__(self, fields):
		self.chrx    = fields[0]
		self.start   = int(fields[3])
		self.stop    = int(fields[4])
		self.exonnum = int(fields[17].replace('"', "")[:-1])
		self.genesym = fields[19].replace('"', "")[:-1]
		self.tr      = fields[13].replace('"', "")[:-1]
		self.count   = 0

def incAllExons(pos, length, array, lowindex, hiindex):
	if hiindex - lowindex < 2:
		# found
		for i in range(lowindex, len(array)):
			if array[i].start < pos + length and array[i].stop > pos:
				array[i].count += 1
			if array[i].start > pos + length:
				break
	else:
		midindex = lowindex + (hiindex - lowindex) / 2
		if array[lowindex].start < pos + length and array[midindex].stop > pos:
			incAllExons(pos, length, array, lowindex, midindex)
		else:
			incAllExons(pos, length, array, midindex, hiindex)

gtf = dict()
gtffile = gzip.open(sys.argv[1])
for line in gtffile:
	if line.startswith("#"):
		continue
	fields = line.split()
	if fields[2] != "exon":
		continue
	exon = Exon(fields)
	if fields[0] not in gtf:
		gtf[fields[0]] = list()
	gtf[fields[0]].append(exon)
gtffile.close()

for i in gtf:
	gtf[i] = sorted(gtf[i], key = lambda x:x.start)

readcount = 0
for i in sys.stdin:
	fields = i.rstrip().split("\t")
	if fields[4] == "255":
		pos    = int(fields[3])
		chrx   = fields[2]
		length = len(fields[9])
		readcount += 1
		incAllExons(pos, length, gtf[chrx], 0, len(gtf[chrx]) - 1)

for i in gtf:
	for exon in gtf[i]:
		print "%s_%s:%d\t%d" % (exon.genesym, exon.tr, exon.exonnum, exon.count)
print "__readcount\t" + str(readcount)
