#!/usr/bin/python3
# 
# given a -m8CB file with exon annotations for the query and subject,
# adjust the subject exon names to match the query exon names
#
# see test_py.sh for sample use
#
################################################################
# copyright (c) 2018 by William R. Pearson and The Rector &
# Visitors of the University of Virginia */
################################################################
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under this License is distributed on an "AS
# IS" BASIS, WITHOUT WRRANTIES OR CONDITIONS OF ANY KIND, either
# express or implied.  See the License for the specific language
# governing permissions and limitations under the License. 
################################################################

import fileinput
import sys
import re
import argparse
import copy

################
# "domain" class that describes a domain/exon alignment annotation
#
class DomAlign:
    def __init__(self, name, info, color, qstart, qend, sstart, send, raw_score, bit_score, ident, qscore, RXRState, fulltext):
        self.name = name
        self.info = info
        self.color_type = ''
        if (not re.search(r'^\d+$',color)):
            m=re.search(r'^(\d+)([a-z]?\w*)$',color)
            if (m):
                (self.color, self.color_type) = m.groups()
                self.color = int(self.color)
        else:
            self.color = int(color)

        self.q_start = qstart
        self.q_end = qend
        self.s_start = sstart
        self.s_end = send
        self.raw_score = raw_score
        self.bit_score = bit_score
        self.percid = ident
        self.q_score = qscore
        self.rxr = RXRState
        self.idnum = 0
        self.overlap_list = []
        self.info_dom = None
        self.text = fulltext
        self.out_str = ''
        self.over_cnt = 0

    def append_overlap(self, overlap_dict):
        self.overlap_list.append(overlap_dict)

    def __str__(self):
#        return "[%d]name: %s : %i-%i : %i-%i I=%.1f Q=%.1f %s" % (self.idnum, self.name, self.q_start, self.q_end, self.s_start, self.s_end, self.percid, self.q_score, self.rxr)
        return "[%d:%s] %i-%i:%i-%i::%s [over:%d]" % (self.idnum, self.rxr, self.q_start, self.q_end, self.s_start, self.s_end, self.name,len(self.overlap_list))

    def print_bar_str(self):  # checking for 'NADA' 
        if (not self.out_str):
            self.out_str = self.text
        return str("|%s"%(self.out_str))

    def make_bar_str(self):  # create original string from values
        bar_str = "|%s:%d-%d:%d-%d:s=%d;b=%.1f;I=%.3f;Q=%.1f;C=%s%s~%d" % (
            self.rxr, self.q_start, self.q_end, self.s_start, self.s_end,
            self.raw_score, self.bit_score, self.percid, self.q_score, self.name, self.info, self.color)

        if (self.color_type):
            bar_str += self.color_type
        return bar_str

################
# "exonInfo" class describes raw (un-aligned) exons with genome coordinates
#
class exonInfo:
    def __init__(self, name, q_target, p_start, p_end, chrom, d_start, d_end, full_text):
        self.name = name
        self.q_target = q_target
        self.p_start = p_start
        self.p_end = p_end
        self.chrom = chrom
        self.d_start = d_start
        self.d_end = d_end
        self.text = full_text
        self.plus_strand = True
        if (d_start > d_end):
            self.plus_strand = False

    def __str__(self):
        rxr_str = "XD"
        if (self.q_target):
            rxr_str="DX"
        return '|%s:%i-%i:%s{%s:%i-%i}' % (rxr_str, self.p_start, self.p_end, self.name, self.chrom, self.d_start, self.d_end)


# Parses domain annotations after split at '|'
#|RX:1-38:3-40:s=37;b=17.0;I=0.289;Q=15.9;C=exon_1~1
#|RX:39-67:41-69:s=78;b=35.8;I=0.483;Q=68.7;C=exon_2~2
#|XR:1-67:3-69:s=115;b=52.8;I=0.373;Q=116.3;C=exon_1~1
#|RX:68-117:72-113:s=14;b=6.4;I=0.385;Q=0.0;C=exon_3~3
#|XR:68-124:70-119:s=-11;b=0.0;I=0.378;Q=0.0;C=exon_2~2
#|XR:125-167:120-165:s=39;b=17.9;I=0.429;Q=18.5;C=exon_3~3
#|RX:118-176:114-175:s=24;b=11.0;I=0.411;Q=1.5;C=exon_4~4
#|RX:177-200:176-198:s=27;b=12.4;I=0.435;Q=4.0;C=exon_5~5
#|XR:168-200:166-198:s=12;b=5.5;I=0.419;Q=0.0;C=exon_4~4
#
def parse_domain(text):
    # takes a domain in string form, turns it into a domain object
    # looks like: RX:5-82:5-82:s=397;b=163.1;I=1.000;Q=453.6;C=C.Thioredoxin~1
    # could also look like: RX:5-82:5-82:s=397;b=163.1;I=1.000;Q=453.6;C=C.Thioredoxin{PF012445}~1

    # get RX/XR and qstart/qstop sstart/sstop as strings
    m = re.search(r'^(\w+):(\d+)-(\d+):(\d+)-(\d+):',text)
    if (m):
        (RXRState, qstart_s, qend_s, sstart_s, send_s) = m.groups()
    else:
        sys.stderr.write("could not parse location: %s\n"%(text))

    # get score, bits, identity, Q info
    m = re.search(r's=(\-?\d+);b=(\-?[\d\.]+);I=([\d\.]+);Q=(\-?\d+\.\d*);',text)
    if (m):
        (r_score_s, b_score_s, ident_s, qscore_s) = m.groups()
    else:
        sys.stderr.write("Error: no scores: %s\n" %(text))
        r_score_s = b_score_s = qscore_s = "-1.0"

    # get domain name/color (and possibly {info})

    (name, color_s) = re.search(r';C=([^~]+)~(.+)$',text).groups()
    info_s=""

    if (re.search(r'\}$',name)):
        (name, info_s) = re.search(r'([^\{]+)(\{[^\}]+\})$',name).groups()

    dom_align = DomAlign(name, info_s, color_s, int(qstart_s), int(qend_s), int(sstart_s), int(send_s), 
                         int(r_score_s), float(b_score_s), float(ident_s),float(qscore_s), RXRState, text)

    return dom_align

# dom_info is like domain, but no scores
################
# exon_info is like domain, but no scores
#
def parse_exon_info(text):
    # takes a domain in string form, turns it into a domain object
    # looks like: DX:1-100;C=C.Thioredoxin~1

    (RXRState, start_s, end_s,name, color) = re.search(r'^(\w+):(\d+)-(\d+);C=([^~]+)~(.*)$',text).groups()
    info = ""
    if (re.search(r'\}$',name)):
        (name, info) = re.search(r'([^\{]+)(\{[^\}]+\})$',name).groups()

    gene_re = re.search(r'^\{([\w\.]+):(\d+)\-(\d+)\}',info)
    if (gene_re):
        (chrom, d_start, d_end) = gene_re.groups()
    else:
        (chrom, d_start, d_end) = ('',-1,-1)
#        sys.stderr.write("genome info not found: %s\n" % (text))

    q_target = True;
    if (RXRState == 'XD'):
        q_target = False

    exon_info = exonInfo(name, q_target, int(start_s), int(end_s), chrom, int(d_start), int(d_end), text)

    return exon_info

def overlap_fract(qdom, sdom):
    # checks if a query and subject domain overlap
    # if they do, return the amount of overlap with respect to each domain
    # how much of query is covered by subject, how much of subject is covered by query

    q_overlap = 0.0
    s_overlap = 0.0

    qq_len = qdom.q_end-qdom.q_start+1   # query alignment length in query coordinates
    qs_len = qdom.s_end-qdom.s_start+1   # query alignment length in subj coordinates
    sq_len = sdom.q_end-sdom.q_start+1   # subj alignment length in query coordinates
    ss_len = sdom.s_end-sdom.s_start+1   # subj alignment length in subject coordinates

    case = -1

    # case (0) no overlap at all
    if (qdom.q_end < sdom.q_start or sdom.s_end < qdom.s_start or qdom.q_start > sdom.q_end or sdom.q_start > qdom.q_end) :
        case = 0
        q_overlap = s_overlap = 0.0
    # case (1) query surrounds subject
    elif (qdom.q_start <= sdom.q_start and qdom.q_end >= sdom.q_end):
        case = 1
        s_overlap = 1.0
        q_overlap = float(sq_len)/qq_len
    # case (2) subject surrounds query
    elif (sdom.s_start <= qdom.s_start and sdom.s_end >= qdom.s_end):
        case = 2
        q_overlap = 1.0
        s_overlap = float(qs_len)/ss_len
    # case (3) query left of subject
    elif (qdom.q_start <= sdom.q_start and qdom.q_end <= sdom.q_end):
        case = 3
        q_overlap = float(qdom.q_end-sdom.q_start+1)/qq_len
        s_overlap = float(qdom.s_end-sdom.s_start+1)/ss_len
    # case (4) subject of left of query
    elif (sdom.s_start <= qdom.s_start and sdom.s_end <= qdom.s_end):
        case = 4
        q_overlap = float(sdom.q_end-qdom.q_start+1)/qq_len
        s_overlap = float(sdom.s_end-qdom.s_start+1)/ss_len

    if (q_overlap > 1.0 or s_overlap > 1.0):
        if (1):
            sys.stderr.write("***%i: qdom: %s sdom: %s\n"% (case,str(qdom),str(sdom)))
            sys.stderr.write(" ** qover %.3f sover: %.3f\n"% (q_overlap, s_overlap))
            sys.stderr.write(" ** qq_len: %d qs_len: %d ss_len: %d sq_len %d\n"%(qq_len, qs_len, ss_len, sq_len))

    return (q_overlap, s_overlap)

####
# parse_protein(result_line)
# takes a protein in string format, turns it into a dictionary properly
# looks like:   sp|P30711|GSTT1_HUMAN   up|Q2NL00|GSTT1_BOVIN   86.67   240     32      0       1       240     1       240     1.4e-123        444.0   16VI7DR6IT3IR15KQ3AI6TI11TA7YH8RC12TA3SN10FL10QETM2AT6VMTA2LV2DG4ND6PS24EK6TA11DV14FSPQ5IL3LMML1WK5RQ   |XR:4-76:4-76:s=327;b=134.6;I=0.895;Q=367.8;C=C.Thioredoxin~1|RX:5-82:5-82:s=356;b=146.5;I=0.902;Q=403.3;C=C.Thioredoxin~1|RX:83-93:83-93:s=52;b=21.4;I=0.818;Q=30.9;C=NODOM~0|XR:77-93:77-93:s=86;b=35.4;I=0.882;Q=72.6;C=NODOM~0|RX:94-110:94-110:s=88;b=36.2;I=0.882;Q=75.0;C=vC.GST_C~2v|XR:94-110:94-110:s=88;b=36.2;I=0.882;Q=75.0;C=vC.GST_C~2v|RX:111-201:111-201:s=409;b=168.3;I=0.868;Q=468.3;C=C.GST_C~2|XR:111-201:111-201:s=409;b=168.3;I=0.868;Q=468.3;C=C.GST_C~2|RX:202-240:202-240:s=154;b=63.4;I=0.795;Q=155.9;C=NODOM~0|XR:202-240:202-240:s=154;b=63.4;I=0.795;Q=155.9;C=NODOM~0
#
# returns [data[x] for x in fields] but also data['q/s_dom_list'] and data['q/sinfo_list']
def parse_protein(line_data,fields, req_name):
    # last part (domain annotions) split('|') and parsed by parse_domain()

    data = {}
    data = dict(list(zip(fields, line_data)))
    if (re.search(r'\|',data['qseqid'])):
        data['qseq_acc'] = data['qseqid'].split('|')[1]
    else:
        data['qseq_acc'] = data['qseqid']

    if (re.search(r'\|',data['sseqid'])):
        data['sseq_acc'] = data['sseqid'].split('|')[1]
    else:
        data['sseq_acc'] = data['sseqid']

    Qdom_list = []
    Sdom_list = []

    Qinfo_list = []
    Sinfo_list = []

    counter = 0

    if ('dom_annot' in data and len(data['dom_annot']) > 0):
        for dom_str in data['dom_annot'].split('|')[1:]:
            if (req_name and not re.search(req_name, dom_str)):
                continue

            counter += 1
            dom = parse_domain(dom_str)
            dom.idnum = counter
            if (dom.rxr == 'RX'):
                Qdom_list.append(dom)
            else:
                Sdom_list.append(dom)

    data['qdom_list'] = Qdom_list
    data['sdom_list'] = Sdom_list

    if ('dom_info' in data and len(data['dom_info']) > 0):
        for info_str in data['dom_info'].split('|')[1:]:
            if (req_name and not re.search(req_name, info_str)):
                continue
            if (not re.search(r'^[DX][XD]',info_str)):
                continue

            dinfo = parse_exon_info(info_str)

            if (dinfo.q_target):
                Qinfo_list.append(dinfo)
            else:
                Sinfo_list.append(dinfo)


        # put links to info_list into dom_list so info_list names can
        # be changed  -- give S/Qinfo's the S/Qdom ids of the overlapping domain

        find_info_overlaps(Qinfo_list, Qdom_list)
        find_info_overlaps(Sinfo_list, Sdom_list)
        
    data['qinfo_list'] = Qinfo_list
    data['sinfo_list'] = Sinfo_list

    return data

################
# "domain" : RX:1-38:3-40:s=37;b=17.0;I=0.289;Q=15.9;C=exon_1~1
# "name"   : like exon_2
# expanded for domain: RX:1-38:3-40:s=37;b=17.0;I=0.289;Q=15.9;C=exon_1{chr1:12345678-123456987}~1
#
def replace_name(domain_text, new_name, new_color_s):
    out = "=".join(domain_text.split("=")[:-1])  # out has everything to last '='

    old_name = domain_text.split(";C=")[-1]
    old_info=""
    
    if (re.search(r'\}~',old_name)):
        (old_info)=re.search(r'(\{[^\}]+\})~',old_name).group(1)

    if (not re.match(r'\d+',new_color_s)):
        new_color_s="0"
    out += "="+new_name+old_info+"~"+new_color_s                 # put it together
    return out

################
# check for overlaps using mid-point
#
def mid_overlaps(qdom_list, sdom_list):

    if (len(qdom_list) != len(sdom_list)):
        return False

    for ix, q_dom in enumerate(qdom_list):
        s_dom = sdom_list[ix]
        q_mid = q_dom.q_start + (q_dom.q_end - q_dom.q_start + 1)/2.0
        if not (q_mid >= s_dom.q_start and q_mid <= s_dom.q_end):
            return False

        q_qfract, q_sfract = overlap_fract(q_dom, s_dom)  # overlap from query perspective
        s_sfract, s_qfract = overlap_fract(s_dom, q_dom)  # overlap from subject perspective

        q_dom.overlap_list.append({"dom": s_dom, "q_over": q_qfract, "s_over": q_sfract})
        s_dom.overlap_list.append({"dom": q_dom, "q_over": s_qfract, "s_over": s_sfract})

    return True

################
# find_overlaps -- populates dom.overlap_list for qdoms, sdoms
#
def find_overlaps(qdom_list, sdom_list, over_thresh):
    # find qdom, sdom overlaps in O(N) time
    #

    if (len(sdom_list) == 0 or len(qdom_list)==0):
        return

    if (len(sdom_list) == len(qdom_list)):  # same number of domains
        if (mid_overlaps(qdom_list, sdom_list)):
            return;
        else:
            for d in qdom_list:
                d.overlap_list = []
            for d in sdom_list:
                d.overlap_list = []

        
    qdom_queue = [x for x in qdom_list]	# build a duplicate list
    sdom_queue = [x for x in sdom_list]

    qdom = qdom_queue.pop(0)	# get the first element of each
    sdom = sdom_queue.pop(0)

    while (True):
        pop_s = pop_q = False

        q_qfract, q_sfract = overlap_fract(qdom, sdom)  # overlap from query perspective
        if (q_qfract > over_thresh  or q_sfract > over_thresh):
            qdom.append_overlap({"dom": sdom, "q_over": q_qfract, "s_over": q_sfract})
            qdom.over_cnt += 1

        s_sfract, s_qfract = overlap_fract(sdom, qdom)  # overlap from query perspective
        if (s_qfract > over_thresh  or s_sfract > over_thresh):
            sdom.append_overlap({"dom": qdom, "q_over": s_qfract, "s_over": s_sfract})
            sdom.over_cnt += 1

        # check to see if we've used up the domain
        if (qdom.s_end >= sdom.s_end):
            pop_s = True
        # else there are more s_dom's that are part of this q_dom

        if (sdom.q_end >= qdom.q_end):
            pop_q = True
        # else there are more q_dom's that are part of this s_dom

#        print 'QS: %s %s\t%s %s' %(pop_q, pop_s, qdom, sdom)

        if (len(qdom_queue) > 0):
            if (pop_q):	# done with this qdom, get next
                qdom = qdom_queue.pop(0)
        elif (pop_q):  # don't break until we try to get the next domain
            break;

        if (len(sdom_queue) > 0):
            if (pop_s):	# done with this sdom, get next
                sdom = sdom_queue.pop(0)
        elif (pop_s):  # don't break until we try to get the next domain
            break;
    ####
    # all done with overlaps

    # # print "overlaps done"
    # for qd in qdom_list:
    #     print qd, qd.over_cnt
    #     for sd in qd.overlap_list:
    #         print " s: q_over %.3f s_over: %.3f %s" % (sd['q_over'], sd['s_over'], str(sd['dom']))
    # print "===="

    # for sd in sdom_list:
    #     print sd, sd.over_cnt
    #     for qd in sd.overlap_list:
    #         print " q: q_over %.3f s_over: %.3f %s" % (qd['q_over'], qd['s_over'], str(qd['dom']))
    # print "===="

################
# info_overlaps -- populates dom.overlap_list for qdoms, sdoms
#
def find_info_overlaps(info_list, dom_list):

    if (len(info_list) == 0 or len(dom_list)==0):
        return

    info_queue = [x for x in info_list]	# build a duplicate list
    dom_queue = [x for x in dom_list]

    info = info_queue.pop(0)	# get the first element of each
    dom = dom_queue.pop(0)

    while (True):
        pop_d = pop_i = False

        if (dom.rxr == 'RX'):  # use dom.q_start/q_end
            if (dom.q_end < info.p_start):
                pop_d = True
            elif (dom.q_end >= info.p_start and dom.q_start <= info.p_end):  # overlap
                dom.info_dom = info
                pop_d = True
                pop_i = True
            elif (info.p_end < dom.q_start):
                pop_i = True
            
        else:			# use dom.s_start/s_end
            if (dom.s_end < info.p_start):
                pop_d = True
            elif (dom.s_end >= info.p_start and dom.s_start <= info.p_end):  # overlap
                dom.info_dom = info
                pop_d = True
                pop_i = True
            elif (info.p_end < dom.s_start):
                pop_i = True

        if (len(info_queue) > 0):
            if (pop_i):	# done with this info, get next
                info = info_queue.pop(0)
        elif (pop_i):  # don't break until we try to get the next domain
            break;

        if (len(dom_queue) > 0):
            if (pop_d):	# done with this dom, get next
                dom = dom_queue.pop(0)
        elif (pop_d):
            break;

################
# build_multi_dict -- builds of dictionaries of multiple overlaps in
#                     qdom.overlap_list or sdom.overlap_list
# returns multi_dict[idnum]
#
def build_multi_dict(dom_list):
    # this code looks for xdom's that are associated with multiple ydoms
    #
    multi_dict = {}  # dict of {qids:/sdom:/qdoms:[]}
    for dom in dom_list:   # for each subject domain
        if (dom.over_cnt <= 1):
            continue

        multi_id_list = []
        multi_dom_list = []
        multi_q_cnt = 0
        for xd_over_yd in dom.overlap_list:  	      # a set of q_doms that overlap the subject
            multi_q_cnt += 1
            multi_id_list.append(xd_over_yd["dom"].idnum)  # these are q_dom idnum's
            multi_dom_list.append(xd_over_yd["dom"])          # these are q_doms

        if (multi_q_cnt > 1):    # only save when two (or more) overlaps
            multi_dict[dom.idnum] = {"yids": multi_id_list, "ydoms":multi_dom_list, 'xdom':dom}

    # # print out current multi_q_list
    # print "--- multi_q dict ---"
    # for db in multi_dict.keys():
    #     print "sdom: %s"%(db)
    #     for ix, qd in enumerate(multi_dict[db]['ydoms']):
    #         print " %d %d: %s"%(ix, multi_dict[db]['yids'][ix], qd)

    # print "--- multi_dict done"

    return multi_dict

################
# find_best_id() -- returns id of domain with longest 'q_over'
#
def find_best_id(overlap_list, over_type):

    max_fract = 0.0
    max_idnum = 0
    for over_d in overlap_list:
        if (over_d[over_type] > max_fract):
            max_idnum = over_d['dom'].idnum
            max_fract = over_d[over_type]

    return max_idnum

################################################################
# final labeling routine -- leave qdom's alone, modify sdoms based on qdoms.
################
# sdom's in more than one qdom are in multi_q_dict[]
# qdom's in more than one sdom are in multi_s_dict[]
# everyone else just gets the qdom name
# returns sdom_displayed_dict{idnum} -- the set of sdoms that have been modified
#
# 13-Nov-2018 -- ensure that there is an info_dom before replacing info_dom.text
#
def label_doms(qdom_list, sdom_list, multi_q_dict, multi_s_dict):

    sdom_displayed_dict = {}
    for qdom in qdom_list:
        # qdom's stay the same
        qdom.out_str = qdom.text

        # check for s_doms with multiple q_doms
        if (qdom.idnum in multi_s_dict):
            # find the best, name it exon_X, find the rest, name it qdom.name
            multi_s_entry = multi_s_dict[qdom.idnum]
            best_id = find_best_id(qdom.overlap_list,'q_over')  # find sdom with most overlap
            for s_over in qdom.overlap_list:   # find the sdom's that overlap this qdom
                sdom = s_over['dom']
                if (sdom.idnum == best_id):
                    sdom.out_str = replace_name(sdom.text, qdom.name, str(qdom.color))
                    if (sdom.info_dom):
                        sdom.info_dom.out_str = replace_name(sdom.info_dom.text,qdom.name, str(qdom.color))
                else:
                    sdom.out_str = replace_name(sdom.text, "exon_X","0")
                    if (sdom.info_dom):
                        sdom.info_dom.out_str = replace_name(sdom.info_dom.text,"exon_X","0")
                sdom_displayed_dict[sdom.idnum] = sdom;
            continue	# prevents re-labeling later

        # check for q_doms with multiple doms
        for sd_over in qdom.overlap_list:
            sdom = sd_over['dom']
            # it might make sense to do this in a second for loop after
            # all the multiple stuff is done
            if (sdom.idnum not in multi_q_dict):
                # this is the simplest case -- sdom.text gets qdom.name
                if (sdom.idnum not in sdom_displayed_dict):
                    sdom.out_str = replace_name(sdom.text, qdom.name, str(qdom.color))
                    if (sdom.info_dom):
                        sdom.info_dom.out_str = replace_name(sdom.info_dom.text,qdom.name, str(qdom.color))
            else:
                # this sdom belongs to multiple q_doms, add each of those q_doms to the name
                exon_str='exon_'
                # "ydoms" here are the qdoms overlapped by sdom
                exon_str += '/'.join([ x.name.split("_")[1] for x in multi_q_dict[sdom.idnum]['ydoms']])
                sdom.out_str = replace_name(sdom.text, exon_str,"0")
                if (sdom.info_dom):
                    sdom.info_dom.out_str = replace_name(sdom.info_dom.text,exon_str,"0")

            sdom_displayed_dict[sdom.idnum]=sdom

    # done with labeling sdoms based on qdoms, but some may be unlabeled
    # check for missing s_doms
    while (len(list(sdom_displayed_dict.keys())) < len(sdom_list)):
        for sdom in sdom_list:
            if (sdom.idnum not in sdom_displayed_dict):
                sdom.out_str = replace_name(sdom.text, "exon_X","0")
                if (sdom.info_dom):
                    sdom.info_dom.out_str = replace_name(sdom.info_dom.text,"exon_X","0")

                sdom_displayed_dict[sdom.idnum] = sdom

    return sdom_displayed_dict

################
#
# aa_to_exon()  --- given a coordinate and the corresponding exon map, return the exon coordinate
# (can only be done for aligned exons)
#
# this version of the function must use an info_list, not an
# align_list, because it uses p_start/p_end rather than q_start/s_start, etc.
# a version using qp_start/sp_start would also need a target argument
#
def aa_to_exon(aa_coords, exon_info_list):

    sorted_aa_coords = sorted(aa_coords)

    pos_strand = True
    if (exon_info_list[0].d_start > exon_info_list[0].d_end):
        pos_strand = False

    ex_x = 0
    exon_coords = []

    aap_x = 0
    this_aap = sorted_aa_coords[aap_x]
    while (ex_x < len(exon_info_list)):
        this_exon = exon_info_list[ex_x]
        if (this_aap <= this_exon.p_end and this_aap >= this_exon.p_start):
            aa_dna_offset = (this_aap - this_exon.p_start) * 3 

            if (pos_strand):
                aa_dna_pos = this_exon.d_start + aa_dna_offset
            else:
                aa_dna_pos = this_exon.d_start - aa_dna_offset

            exon_coords.append({'chrom':this_exon.chrom, 'dpos':aa_dna_pos})
            aap_x += 1
            if (aap_x < len(sorted_aa_coords)):
                this_aap = sorted_aa_coords[aap_x]
            else:
                break
        else:
            ex_x += 1

    aa_coord_dict = {}
    for aap_x, aap in enumerate(sorted_aa_coords):
        aa_coord_dict[aap] = exon_coords[aap_x]

    return [aa_coord_dict[ax] for ax in aa_coords]

################
# 
def set_data_fields(args, line_data) :

    field_str = 'qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore BTOP dom_annot'
    field_qs_str = 'qseqid qlen sseqid slen pident length mismatch gapopen qstart qend sstart send evalue bitscore BTOP dom_annot'

    if (len(line_data) > 1) :
        if ((not args.have_qslen) and  re.search(r'\d+',line_data[1])):
            args.have_qslen=True

        if ((not args.dom_info) and re.search(r'^\|[DX][XD]\:',line_data[-1])):
            args.dom_info = True

    end_field = -1
    fields = field_str.split(' ')

    if (args.have_qslen):
        fields = field_qs_str.split(' ')

    if (args.dom_info):
        fields.append('dom_info')
        end_field = -2

    return (fields, end_field)

################################################################
#
# main program 
# print "#"," ".join(sys.argv)

def main():

    parser=argparse.ArgumentParser(description='rename_exons.py result_file.m8CB : re-label subject exons to match query')
    parser.add_argument('--have_qslen', help='bl_tab fields include query/subject lengths',dest='have_qslen',action='store_true',default=False)
    parser.add_argument('--dom_info', help='raw domain coordinates included',action='store_true',default=False)
    parser.add_argument('--fill_gcoords', help='fill in genomic coordinates',action='store_true',default=False)
    parser.add_argument('files', metavar='FILE', nargs='*', help='files to read, if empty, stdin is used')

    args=parser.parse_args()

    end_field = -1
    data_fields_reset=False

    (fields, end_field) = set_data_fields(args, [])

    if (args.have_qslen and args.dom_info):
        data_fields_reset=True

    saved_qdom_list = []
    qdom_list = []

    for line in fileinput.input(args.files):
    # pass through comments
        if (line[0] == '#'):
            print(line, end='')	# ',' because have not stripped
            continue

        ################
        # break up tab fields, check for extra fields
        line = line.strip('\n')
        line_data = line.split('\t')
        if (not data_fields_reset):     # look for --have_qslen number, --dom_info data, even if not set
            (fields, end_field) = set_data_fields(args, line_data)
            data_fields_reset = True

        ################
        # get exon annotations
        data = parse_protein(line_data,fields,"exon")	# get score/alignment/domain data

        if (len(data['sdom_list'])==0 and len(data['qdom_list'])==0):
            print(line)	# no domains to be edited, print stripped line and contine
            continue

        # qdom_list=[] outside of loop for cases where the qseqid==sseqid match is not first
        if len(data['qdom_list'])== 0:
            if data['qseqid'] == data['sseqid']:
                saved_qdom_list = [ copy.deepcopy(x) for x in data['sdom_list']]
                max_sdom_id=len(data['sdom_list'])+1
                for qdom in saved_qdom_list:
                    qdom.rxr = 'RX'
                    qdom.idnum = max_sdom_id
                    max_sdom_id += 1

            qdom_list = [copy.deepcopy(x) for x in saved_qdom_list]
        else:
            qdom_list = data['qdom_list']

        # print out non-exon info
    
        if (len(qdom_list) == 0):
            print(line)
            continue

        btab_str = '\t'.join(str(data[x]) for x in fields[:end_field])
        # print  # comment out for single line

        ################
        # find overlaps and multi-overlaps
        #
        find_overlaps(qdom_list,data['sdom_list'], 0.2)

        multi_q_dict = build_multi_dict(data['sdom_list'])  # keys are sdoms hitting multiple qdoms
        multi_s_dict = build_multi_dict(qdom_list)  # keys are qdoms hitting mulitple sdoms

        ################
        # label qdoms, relabel sdoms
        #
        sdom_displayed_dict = label_doms(qdom_list, data['sdom_list'], multi_q_dict, multi_s_dict)

        ################
        # print exon annotations
        #
        q_exon_list = data['qdom_list']

        s_exon_list = [sdom_displayed_dict[x] for x in list(sdom_displayed_dict.keys())]

        ################
        # if args.fill_gcoords, then do the transformations on the current exon lists
        
        if (args.fill_gcoords):
            sa_from_qa = []
            for q_ex in q_exon_list:
                sa_from_qa.append(q_ex.q_start)
                sa_from_qa.append(q_ex.q_end)

            # have list of coordinates, map them to exon
            sex_from_qa2sa = aa_to_exon(sa_from_qa,data['sinfo_list'])

            for iqx, q_ex in enumerate(q_exon_list):
                sg_start = sex_from_qa2sa[2*iqx]
                sg_end = sex_from_qa2sa[2*iqx+1]
                sg_replace="::%s:%d-%d}"%(sg_start['chrom'],sg_start['dpos'],sg_end['dpos'])
                q_ex.text=re.sub(r'\}',sg_replace,q_ex.text)
                q_ex.out_str=re.sub(r'\}',sg_replace,q_ex.out_str)
                
            qa_from_sa = []
            for s_ex in s_exon_list:
                qa_from_sa.append(s_ex.q_start)
                qa_from_sa.append(s_ex.q_end)

            # have list of coordinates, map them to exon
            qex_from_sa2qa = aa_to_exon(qa_from_sa,data['qinfo_list'])

            for isx, s_ex in enumerate(s_exon_list):
                qg_start = sex_from_qa2sa[2*iqx]
                qg_end = sex_from_qa2sa[2*iqx+1]
                qg_replace="{%s:%d-%d::"%(sg_start['chrom'],sg_start['dpos'],sg_end['dpos'])
                s_ex.text=re.sub(r'\{',qg_replace,s_ex.text)
                s_ex.out_str = re.sub(r'\{',qg_replace,s_ex.out_str)

        sorted_exon_list = sorted(q_exon_list+s_exon_list,key = lambda r: r.idnum)

        dom_bar_str = ''
        for exon in sorted_exon_list:
            # print exon.print_bar_str()     # for multi-line output
            dom_bar_str += exon.print_bar_str() 

        info_bar_str = ''
        for info in data['qinfo_list'] + data['sinfo_list']:
            info_bar_str += info.text 

        print('\t'.join((btab_str, dom_bar_str, info_bar_str)))

################
# run the program ...

if __name__ == '__main__':
    main()