##################################
#                                #
# Last modified 01/29/2009       # 
#                                #
# Georgi Marinov                 #
#                                # 
##################################

import sys
import math

try:
	import psyco
	psyco.full()
except:
	pass

def run():

    if len(sys.argv) < 6:
        print 'usage: python %s region [labelFields (0,1,2,3,...)] chr_field peakfield refFlat outfilename [-noPeak]' % sys.argv[0]
        sys.exit(1)

    ERANGE = sys.argv[1]
    labels=sys.argv[2].split(',')
    labelFields=[]
    for i in labels:
        labelFields.append(int(i))
    labelFields.sort()
    ERANGE_chrfieldID=int(sys.argv[3])
    ERANGE_peakfieldID=int(sys.argv[4])
    refFlat = sys.argv[5]
    outputfilename = sys.argv[6]

    noPeak=False
    if '-noPeak' in sys.argv:
        noPeak=True

    outfile = open(outputfilename, 'w')

    input_stream = open(refFlat)
    TSSDict={}
    for line in input_stream:
        fields=line.strip().split('\t')
        chr=fields[2]
        orientation=fields[3]
        left=int(fields[4])
        right=int(fields[5])
        name=fields[0]
        if orientation == '+':
            TSS=left
        if orientation == '-':
            TSS=right
        if TSSDict.has_key(chr):
            TSSDict[chr].append((TSS,orientation,name))
        else:
            TSSDict[chr]=[]
            TSSDict[chr].append((TSS,orientation,name))
            
    input_stream = open(ERANGE)
    for line in input_stream:
        if line[0]=='#':
            continue
        fields=line.strip().split('\t')
        if noPeak:
            left=int(fields[ERANGE_chrfieldID+1])
            right=int(fields[ERANGE_chrfieldID+2])
            peak=int((left+right)/2.0)
        else:
            peak=int(fields[ERANGE_peakfieldID])
        chr=fields[ERANGE_chrfieldID]
        if TSSDict.has_key(chr):
            pass
        else:
            continue
        distanceList=[]
        for (TSS,orientation,name) in TSSDict[chr]:
            AbsoluteDistance=math.fabs(TSS-peak)
            if orientation=='+':
                distance=peak-TSS
            if orientation=='-':
                distance=TSS-peak
            distanceList.append((AbsoluteDistance,name,distance))
        distanceList.sort()
        outline=''
        for ID in labelFields:
            outline=outline+fields[ID]+'\t'
        outline=outline + distanceList[0][1] + '\t' + str(distanceList[0][2]) + '\t' + distanceList[1][1] + '\t' + str(distanceList[1][2])
        outfile.write(outline+'\n')
    
    outfile.close()

run()