<gh_stars>100-1000
#!/usr/bin/python
# Copyright (C) 2010 Ion Torrent Systems, Inc. All Rights Reserved

import sys
import bisect
from collections import defaultdict
from operator import itemgetter

region_ids = defaultdict(lambda: defaultdict(defaultdict))
leftAlignedHotspotStarts = defaultdict(lambda: defaultdict(defaultdict))

def getLociDataNew(chrom, pos):
	if have_bed:
		ids = region_ids[chrom][pos]
		if ids != None and len(ids) > 0:
			return ids
		else:
    			return ("N/A","N/A","N/A")

# assumes a BED file sorted by chrom,start,end
have_bed = (len(sys.argv) >= 5) and (sys.argv[4] != "") and (sys.argv[3] != "")

if have_bed:
	leftAlignedGenes = open(sys.argv[3], 'r')
	for lines in leftAlignedGenes:
		if len(lines) == 0:
			continue
		tmp = lines.split('\t')
		chrom = tmp[0].strip()
        	if chrom.startswith('track '):
            		continue
        	# avoid 0-length BED regions (un-anchored inserts)
        	rsrt = int(tmp[1])+1
        	rend = int(tmp[2])
        	if rsrt > rend:
            		continue
		hotspotId = tmp[3].strip()
		leftAlignedHotspotStarts[hotspotId] = rsrt
	leftAlignedGenes.close()
		
if have_bed:
    genes = open(sys.argv[4],'r')
    for lines in genes:
        if len(lines) == 0:
            continue
        tmp = lines.split('\t')
        chrom = tmp[0].strip()
        if chrom.startswith('track '):
            continue
        # avoid 0-length BED regions (un-anchored inserts)
        rsrt = int(tmp[1])+1
        rend = int(tmp[2])
	hotspotId = 'N/A'
	if len(tmp) > 3 :
		hotspotId = tmp[3].strip()
		if(leftAlignedHotspotStarts[hotspotId] != rsrt):
			hotspotId = hotspotId + "*";
        if rsrt > rend:
            continue
        while (rsrt <= rend):
		if len(tmp) < 6:
            		region_ids[chrom].setdefault(rsrt,[]).append( ('N/A', 0, 'N/A'))
        	else:
            		region_ids[chrom].setdefault(rsrt,[]).append((tmp[len(tmp)-1].strip(),tmp[len(tmp)-2].strip(),hotspotId) )
		rsrt += 1
    genes.close()

# assumes input file has forward+reverse read coverage (Uppercase fields) and reverse coverage (Lowercase fields)
inf = open(sys.argv[1],'r')
out = open(sys.argv[2],'w')
out.write("Chrom\tPosition\tTarget ID\tHotSpot ID\tRef\tCov\tA Reads\tC Reads\tG Reads\tT Reads\t+Cov\t-Cov\tDeletions\n")
for lines in inf:
    if lines[0]=='#':
        continue
    else:
	regIds = []
	hotspotIds = []
        (contig,position,ref,cov,cov_A,cov_C,cov_G,cov_T,cov_D,cov_a,cov_c,cov_g,cov_t,cov_d) = lines.split('\t')
	startPosition = int(position)
	ids = getLociDataNew( contig, startPosition )
	for item in ids:
		if(not item[0] in regIds):
			regIds.append(item[0])
		if(not item[2] in hotspotIds):
			hotspotIds.append(item[2])
	#(region_id,aux_id,hotspot_id) = getLociData( contig, int(position) )
        # forward/reverse coverage does not include coverage by deletions
        cov_r = int(cov_a) + int(cov_c) + int(cov_g) + int(cov_t)
        cov_f = int(cov) - int(cov_D) - cov_r
        
	out.write("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%d\t%d\t%s\n" % (contig,position, ";".join(regIds), ";".join(hotspotIds),ref,cov,cov_A,cov_C,cov_G,cov_T,cov_f,cov_r,cov_D))
out.close()
inf.close()
