##################################
#                                #
# Last modified 2016/12/30       # 
#                                #
# Georgi Marinov                 #
#                                # 
##################################

import sys
import string
import math

def run():

    if len(sys.argv) < 5:
        print 'usage: python %s regions chrFieldID phastConsWig minScore outputfilename' % sys.argv[0]
        print '\t use - for stdIn for the wig file'
        sys.exit(1)
    
    regions = sys.argv[1]
    chrFieldID = int(sys.argv[2])
    input = sys.argv[3]
    minScore = float(sys.argv[4])
    outfilename = sys.argv[5]

    RegionsCoverageDict = {}
    linelist = open(regions)
    for line in linelist:
        if line.startswith('#'):
            continue
        fields = line.strip().split('\t')
        chr = fields[chrFieldID]
        left = int(fields[chrFieldID + 1])
        right = int(fields[chrFieldID + 2])
        if RegionsCoverageDict.has_key(chr):
            pass
        else:
            RegionsCoverageDict[chr] = {}
        for i in range(left,right):
            RegionsCoverageDict[chr][i] = 0

    if input == '-':
        lineslist  = sys.stdin
    else:
        lineslist = open(input)
    InRegion = False
    for line in lineslist:
        if line.startswith('fixedStep'):
            chr=line.split(' ')[1].split('=')[1]
            currentPos=int(line.split(' ')[2].split('=')[1])
            continue
        score = float(line.split('\n')[0])
        if RegionsCoverageDict.has_key(chr) and RegionsCoverageDict[chr].has_key(currentPos):
            RegionsCoverageDict[chr][currentPos] = score
        currentPos+=1

    outfile = open(outfilename, 'w')

    linelist = open(regions)
    for line in linelist:
        if line.startswith('#'):
            outline = line.strip() + '\tconservation_score\tfraction_conserved'
            outfile.write(outline + '\n')
            continue
        fields = line.strip().split('\t')
        chr = fields[chrFieldID]
        left = int(fields[chrFieldID + 1])
        right = int(fields[chrFieldID + 2])
        consScore = 0.0
        consBPs = 0.0
        for i in range(left,right):
            if RegionsCoverageDict[chr][i] >= minScore:
                consScore += RegionsCoverageDict[chr][i]
                if RegionsCoverageDict[chr][i] > 0:
                    consBPs += 1
        fractionCons = consBPs/(right-left)
        outline = line.strip() + '\t' + str(consScore) + '\t' + str(fractionCons)
        outfile.write(outline + '\n')

    outfile.close()
   
run()