dnarrange

#! /usr/bin/env python
# Author: Martin C. Frith 2018
# SPDX-License-Identifier: GPL-3.0-or-later

from __future__ import print_function

import collections
import functools
import gzip
import heapq
import logging
import math
import optparse
import os
import re
import signal
import sys
import textwrap
from itertools import chain, groupby, islice
from operator import itemgetter

def myOpen(fileName):  # faster than fileinput
    if fileName == "-":
        return sys.stdin
    if fileName.endswith(".gz"):
        return gzip.open(fileName, "rt")  # xxx dubious for Python2
    return open(fileName)

def connectedComponent(adjacencyList, nodePriorities, isNew, i, isFlipped):
    newItem = nodePriorities[i], i, isFlipped
    heap = [newItem]
    isNew[i] = False
    while heap:
        _, j, isFlipped = heapq.heappop(heap)
        yield j, isFlipped
        for k, isOpposite in adjacencyList[j]:
            if isNew[k]:
                newItem = nodePriorities[k], k, isOpposite != isFlipped
                heapq.heappush(heap, newItem)
                isNew[k] = False

def connectedComponents(adjacencyList, nodePriorities, isRev):
    isNew = [True for i in adjacencyList]
    s = sorted((x, i) for i, x in enumerate(nodePriorities))
    for _, i in s:
        if isNew[i]:
            yield list(connectedComponent(adjacencyList, nodePriorities,
                                          isNew, i, isRev[i]))

def adjacencyListFromLinks(numOfNodes, links):
    adjacencyList = [[] for i in range(numOfNodes)]
    for i, j, isOpposite in links:
        adjacencyList[i].append((j, isOpposite))
        adjacencyList[j].append((i, isOpposite))
    return adjacencyList

def dataFromMafFields(fields):
    seqName, start, span, strand, seqLen, seq = fields[1:7]
    beg = int(start)
    seqLen = int(seqLen)
    if strand == "-":
        beg -= seqLen  # use negative coordinates for reverse strands
    return seqName, seqLen, beg, seq

def splitAtBigGaps(opts, refBeg, qryBeg, refSeq, qrySeq):
    minGap = int(math.ceil(opts.min_gap))
    alnLen = len(qrySeq)
    gapString = "-" * min(minGap, alnLen + 1)
    alnPos = 0
    while alnPos < alnLen:
        gapBeg = qrySeq.find(gapString, alnPos)
        if gapBeg < 0:
            gapBeg = gapEnd = alnLen
        else:
            gapEnd = gapBeg + minGap
            while gapEnd < alnLen and qrySeq[gapEnd] == "-":
                gapEnd += 1
            while gapEnd < alnLen and refSeq[gapEnd] == "-":
                gapEnd += 1
            while gapBeg > alnPos and refSeq[gapBeg - 1] == "-":
                gapBeg -= 1
        dist = gapBeg - alnPos
        refEnd = refBeg + dist - refSeq.count("-", alnPos, gapBeg)
        qryEnd = qryBeg + dist - qrySeq.count("-", alnPos, gapBeg)
        yield qryBeg, qryEnd, refBeg, refEnd
        gapLen = gapEnd - gapBeg
        refBeg = refEnd + gapLen - refSeq.count("-", gapBeg, gapEnd)
        qryBeg = qryEnd + gapLen - qrySeq.count("-", gapBeg, gapEnd)
        alnPos = gapEnd

def splitTabAtBigGaps(opts, refBeg, qryBeg, gapText):
    refEnd = refBeg
    qryEnd = qryBeg
    for i in gapText.split(","):
        if ":" in i:
            refInsert, qryInsert = i.split(":")
            refInsert = int(refInsert)
            qryInsert = int(qryInsert)
            if refInsert >= opts.min_gap:
                yield qryBeg, qryEnd, refBeg, refEnd
                refBeg = refEnd + refInsert
                qryBeg = qryEnd + qryInsert
            refEnd += refInsert
            qryEnd += qryInsert
        else:
            size = int(i)
            refEnd += size
            qryEnd += size
    yield qryBeg, qryEnd, refBeg, refEnd

def mismapFromFields(fields):
    for i in fields:
        if i.startswith("mismap="):
            return float(i[7:])
    return 0.0

def alignmentsFromLines(opts, lines):
    refSeq = qrySeq = None
    qryId = 0
    qryEnd = 0
    for line in lines:
        if line[0].isdigit():
            fields = line.split()
            if len(fields) < 12:  # shrunk format
                qryBeg = qryEnd + int(fields[0])
                qryLen = int(fields[1])
                qryEnd = qryBeg + qryLen
                if len(fields) > 4:
                    refName = fields[4]
                    refBeg = int(fields[2])
                else:
                    refBeg = refEnd + int(fields[2])
                refEnd = refBeg + qryLen + int(fields[3])
                alns = [(qryBeg, qryEnd, refBeg, refEnd)]
                yield str(qryId), 0, refName, alns, fields
            else:  # LAST tabular format
                mismap = mismapFromFields(fields[12:])
                if mismap <= opts.max_mismap:
                    refName, refLen, refBeg, rj = dataFromMafFields(fields)
                    qryName, qryLen, qryBeg, qj = dataFromMafFields(fields[5:])
                    alns = splitTabAtBigGaps(opts, refBeg, qryBeg, qj)
                    yield qryName, qryLen, refName, alns, fields
        elif line[0].isalpha():  # MAF format
            if line[0] == "a":
                if qrySeq:
                    yield qryName, qryLen, refName, alns, mafLines
                mafLines = []
                refSeq = qrySeq = None
                mismap = mismapFromFields(line.split())
            elif line[0] == "s" and mismap <= opts.max_mismap:
                fields = line.split()
                if refSeq is None:
                    refName, refLen, refBeg, refSeq = dataFromMafFields(fields)
                else:
                    qryName, qryLen, qryBeg, qrySeq = dataFromMafFields(fields)
                    alns = splitAtBigGaps(opts, refBeg, qryBeg, refSeq, qrySeq)
            mafLines.append(line)
        else:
            qryId += 1
            qryEnd = 0
    if qrySeq:
        yield qryName, qryLen, refName, alns, mafLines

def isCircularChromosome(name):
    return name in ("chrM", "M")  # xxx ???

def isKnownChromosome(name):
    unknownPrefixes = "chrUn", "Un"  # xxx ???
    return not name.startswith(unknownPrefixes)

def chromosomeFromName(name):
    return name.split("_")[0]  # e.g. chr5_random -> chr5

def isDifferentChromosomes(nameX, nameY):
    return (isKnownChromosome(nameX) and isKnownChromosome(nameY) and
            chromosomeFromName(nameX) != chromosomeFromName(nameY))

def refNameAndStrand(alignment):
    return alignment[3], alignment[4] < 0

def knownChromosomes(alignments):
    for i in alignments:
        refName = i[3]
        if isKnownChromosome(refName):
            yield chromosomeFromName(refName)

def isInterChromosome(alignments):
    """Is any pair of alignments on different chromosomes?"""
    return len(set(knownChromosomes(alignments))) > 1

def isInterStrand(alignments):
    """Is any pair of alignments on opposite strands of the same chromosome?"""
    names = set(i[3] for i in alignments)
    namesAndStrands = set(map(refNameAndStrand, alignments))
    return len(namesAndStrands) > len(names)

def isNonlinear(sortedAlignmentsOfOneQuery, opts):
    """Is any pair of alignments non-colinear on the same strand?"""
    maxCoordinates = {}
    for i in sortedAlignmentsOfOneQuery:
        if isCircularChromosome(i[3]):
            continue
        k = refNameAndStrand(i)
        if k in maxCoordinates:
            m = maxCoordinates[k]
            if m >= i[4] + opts.min_rev:
                return True
            if i[5] > m:
                maxCoordinates[k] = i[5]
        else:
            maxCoordinates[k] = i[5]
    return False

def isBigGap(sortedAlignmentsOfOneQuery, opts):
    """Is any pair of adjacent aligments separated by a big genomic gap?"""
    for j, y in enumerate(sortedAlignmentsOfOneQuery):
        if j:
            x = sortedAlignmentsOfOneQuery[j - 1]
            if refNameAndStrand(x) == refNameAndStrand(y):
                if y[4] - x[5] >= opts.min_gap:
                    return True
    return False

def rearrangementType(opts, alignmentsOfOneQuery):
    if opts.insert:
        seqNames = set(i[3] for i in alignmentsOfOneQuery)
        return "I" if len(seqNames) > 1 and opts.insert in seqNames else None
    if "C" in opts.types and isInterChromosome(alignmentsOfOneQuery):
        return "C"
    if "S" in opts.types and isInterStrand(alignmentsOfOneQuery):
        return "S"
    if "N" in opts.types and isNonlinear(alignmentsOfOneQuery, opts):
        return "N"
    if "G" in opts.types and isBigGap(alignmentsOfOneQuery, opts):
        return "G"
    return None

def qryFwdAlns(alignmentGroup):
    for qryName, qryLen, refName, alns, junk in alignmentGroup:
        for qryBeg, qryEnd, refBeg, refEnd in alns:
            if qryBeg < 0:  # use forward strand of query:
                qryBeg, qryEnd = -qryEnd, -qryBeg
                refBeg, refEnd = -refEnd, -refBeg
            yield qryName, qryBeg, qryEnd, refName, refBeg, refEnd

def newAlnFromOldAln(oldAln, qryNum, alnNum):
    qryName, qryBeg, qryEnd, refName, refBeg, refEnd = oldAln
    if refBeg < 0:  # use forward strand of reference:
        refBeg, refEnd = -refEnd, -refBeg
        qryBeg, qryEnd = -qryEnd, -qryBeg
    isRev = (qryBeg < 0)
    return qryNum, qryBeg, qryEnd, refName, refBeg, refEnd, alnNum, [], isRev

def alignmentsPerRearrangedQuerySequence(opts, fileNames):
    qryNum = 0
    alnNum = 0
    for fileNum, fileName in enumerate(fileNames):
        logging.info("reading {0}...".format(fileName))
        alignments = alignmentsFromLines(opts, myOpen(fileName))
        for key, group in groupby(alignments, itemgetter(0, 1)):
            qryName, qryLen = key
            group = list(group)
            alignmentsOfOneQuery = sorted(qryFwdAlns(group))
            rType = rearrangementType(opts, alignmentsOfOneQuery)
            if rType:
                newAlns = []
                for i in alignmentsOfOneQuery:
                    newAlns.append(newAlnFromOldAln(i, qryNum, alnNum))
                    alnNum += 1
                alignmentTexts = [i[4] for i in group]
                yield newAlns, alignmentTexts, fileNum + 1, qryName, rType
                qryNum += 1

def alignedQueryLength(alignmentsOfOneQuery):
    return sum(i[2] - i[1] for i in alignmentsOfOneQuery)

def addNgOverlaps(okAlignmentsInGenomeOrder, ngAlignmentsInGenomeOrder):
    logging.info("finding overlaps for exclusion...")
    n = len(ngAlignmentsInGenomeOrder)
    i = 0
    for alnA in okAlignmentsInGenomeOrder:
        refNameA, refBegA, refEndA = alnA[3:6]
        overlapsA = alnA[7]
        while i < n and ngAlignmentsInGenomeOrder[i][3] < refNameA:
            i += 1
        j = i
        while j < n:
            alnB = ngAlignmentsInGenomeOrder[j]
            if alnB[3] > refNameA or alnB[4] >= refEndA:
                break
            if alnB[5] > refBegA:
                overlapsA.append(alnB[6])
            else:
                ngAlignmentsInGenomeOrder[j] = ngAlignmentsInGenomeOrder[i]
                i += 1
            j += 1

def addOverlaps(myAlignmentsInGenomeOrder):
    logging.info("finding overlaps...")
    stash = []
    refName = ""
    for alnB in myAlignmentsInGenomeOrder:
        if alnB[3] != refName:
            n = 0
            refName = alnB[3]
        qryNum = alnB[0]
        refBeg = alnB[4]
        j = 0
        for alnA in islice(stash, n):
            if alnA[5] > refBeg:  # overlap in ref
                if alnA[0] < qryNum:
                    alnA[7].append(alnB[6])
                elif alnA[0] > qryNum:
                    alnB[7].append(alnA[6])
                stash[j] = alnA
                j += 1
        if len(stash) == j: stash.append(None)
        stash[j] = alnB
        n = j + 1

def delOverlaps(alignments):
    for i in alignments:
        myList = i[7]
        del myList[:]

def overlapsOfOneQuery(alignments, alignmentsOfOneQuery):
    for alnA in alignmentsOfOneQuery:
        alnNumA = alnA[6]
        for alnNumB in alnA[7]:
            qryNumB = alignments[alnNumB][0]
            yield qryNumB, alnNumA, alnNumB

def isAdjacent(alnX, alnY):
    return alnX[6] + 1 == alnY[6] or alnY[6] + 1 == alnX[6]

def isNonlinearPair(opts, types, alnX, alnY):
    if alnX[3] != alnY[3]:
        if "I" in types and opts.insert in (alnX[3], alnY[3]):
            return True
        return "C" in types and isDifferentChromosomes(alnX[3], alnY[3])
    if alnX[8] is not alnY[8]:
        return "S" in types
    if alnX[1] < alnY[1]:
        gap = alnY[4] - alnX[5]
    else:
        gap = alnX[4] - alnY[5]
    if "N" in types and gap + opts.min_rev <= 0:
        if not isCircularChromosome(alnX[3]):
            return True
    if "G" in types and gap >= opts.min_gap:
        if isAdjacent(alnX, alnY):
            return True
    return False

def alignmentEdges(alnA, alnB, isGetEnds):
    if isGetEnds:
        return -alnA[2], -alnA[5], -alnB[2], -alnB[5]
    else:
        return alnA[1], alnA[4], alnB[1], alnB[4]

def isSharedRearrangement(opts, alnAX, alnAY, alnBX, alnBY):
    # alnAX of query sequence A overlaps alnBX of query sequence B
    # alnAY of query sequence A overlaps alnBY of query sequence B
    # alnAX is upstream of alnAY in query sequence A

    qryAX, refAX, qryBX, refBX = alignmentEdges(alnAX, alnBX, not alnAX[8])
    qryAY, refAY, qryBY, refBY = alignmentEdges(alnAY, alnBY, alnAY[8])

    qryDistanceA = qryAX + qryAY
    qryDistanceB = qryBX + qryBY
    begDiff = refAX - refBX
    endDiff = refBY - refAY
    if abs(qryDistanceB - qryDistanceA + begDiff - endDiff) > opts.max_diff:
        return False

    if alnAX[8] is not alnAY[8] or alnAX[3] != alnAY[3]:
        return True

    gapA = refAX + refAY
    gapB = refBX + refBY
    gapAtoB = refAX + refBY
    gapBtoA = refBX + refAY
    gapMin = min(gapA, gapB)
    gapMax = max(gapA, gapB)

    return (gapMax <= -opts.min_rev and gapMax * 2 <= gapMin
            and gapAtoB < 0 and gapBtoA < 0
            or
            gapMin >= opts.min_gap and gapMin * 2 >= gapMax
            and isAdjacent(alnBX, alnBY)
            and gapAtoB > 0 and gapBtoA > 0)

def sharedRearrangement(opts, types, alignmentsA, alignmentsB,
                        overlapsBetweenTwoQueries):
    # "A" refers to a query sequence
    # "B" refers to a different query sequence
    groups = groupby(overlapsBetweenTwoQueries, itemgetter(1))
    overlapsPerAlnA = [(alnNumA, [i[2] for i in v]) for alnNumA, v in groups]

    for alnNumAY, alnNumsBY in overlapsPerAlnA:
        for alnNumAX, alnNumsBX in overlapsPerAlnA:
            if alnNumAX == alnNumAY:
                break
            alnAX = alignmentsA[alnNumAX]
            alnAY = alignmentsA[alnNumAY]
            if not isNonlinearPair(opts, types, alnAX, alnAY):
                continue
            isRevAX = alnAX[8]
            isRevAY = alnAY[8]
            for alnNumBY in alnNumsBY:
                alnBY = alignmentsB[alnNumBY]
                for alnNumBX in alnNumsBX:
                    alnBX = alignmentsB[alnNumBX]
                    if isRevAY is alnBY[8]:
                        if (alnNumBX < alnNumBY and isRevAX is alnBX[8] and
                            isSharedRearrangement(opts,
                                                  alnAX, alnAY, alnBX, alnBY)):
                            return "+"
                    else:
                        if (alnNumBX > alnNumBY and isRevAX is not alnBX[8] and
                            isSharedRearrangement(opts,
                                                  alnAX, alnAY, alnBX, alnBY)):
                            return "-"
    return None

def isNoSharedRearrangement(opts, okAlignments, ngAlignments,
                            alignmentsOfOneOkQuery, rType):
    types = rType if opts.filter > 0 else opts.types
    overlaps = sorted(overlapsOfOneQuery(ngAlignments, alignmentsOfOneOkQuery))
    for qryNumB, g in groupby(overlaps, itemgetter(0)):
        if sharedRearrangement(opts, types, okAlignments, ngAlignments, g):
            return False
    return True

def linksBetweenQueries(opts, alignments, alignmentsPerQuery):
    logging.info("linking...")
    types = "I" if opts.insert else opts.types
    for alignmentsOfOneQuery in alignmentsPerQuery:
        qryNumA = alignmentsOfOneQuery[0][0]
        overlaps = sorted(overlapsOfOneQuery(alignments, alignmentsOfOneQuery))
        for qryNumB, g in groupby(overlaps, itemgetter(0)):
            strand = sharedRearrangement(opts, types,
                                         alignments, alignments, g)
            if strand:
                yield qryNumA, qryNumB, strand == "-"

def linksBetweenClumps(alignments, alignmentsPerQuery, clumps):
    clumpInfoPerQuery = [len(clumps)] * len(alignmentsPerQuery)
    for clumpNum, clump in enumerate(clumps):
        for qryNum, isFlipped in clump:
            clumpInfoPerQuery[qryNum] = clumpNum
    for a, clumpA in enumerate(clumps):
        for qryNumA, isFlippedA in clumpA:
            for alnA in alignmentsPerQuery[qryNumA]:
                for alnNumB in alnA[7]:
                    alnB = alignments[alnNumB]
                    qryNumB = alnB[0]
                    b = clumpInfoPerQuery[qryNumB]
                    if b < len(clumps) and b != a:
                        yield min(a, b), max(a, b), False

def addSharedJumps(opts, jumpsInGenomeOrder):
    stash = []
    n = 0
    for jumpB in jumpsInGenomeOrder:
        alnBX, alnBY, overlapsB = jumpB
        qryNumB = alnBX[0]
        refBegBX = alnBX[4]
        refBegBY = alnBY[4]
        refEndBY = alnBY[5]
        isRevB = (alnBX[6] > alnBY[6])
        if isRevB: alnBX, alnBY = alnBY, alnBX
        j = 0
        for jumpA in islice(stash, n):
            alnAX, alnAY, overlapsA = jumpA
            if alnAX[5] > refBegBX:
                qryNumA = alnAX[0]
                if (alnAY[4] < refEndBY and alnAY[5] > refBegBY and
                    qryNumA != qryNumB):
                    if isRevB: alnAX, alnAY = alnAY, alnAX
                    if isSharedRearrangement(opts, alnBX, alnBY, alnAX, alnAY):
                        overlapsA.append(qryNumB)
                        overlapsB.append(qryNumA)
                stash[j] = jumpA
                j += 1
            else:
                overlapsA[:] = list(set(overlapsA))  # try to save memory
        if len(stash) == j: stash.append(None)
        stash[j] = jumpB
        n = j + 1

def isAllJumpsSupported(opts, goodQryNums, alignmentsOfOneQuery):
    for j, y in enumerate(alignmentsOfOneQuery):
        if j:
            x = alignmentsOfOneQuery[j - 1]
            if isNonlinearPair(opts, opts.types, x, y):
                if len(set(y[7]) & goodQryNums) < opts.min_cov:
                    return False
    return True

def querySortKey(alignmentsOfOneQuery):
    return min(a[3:6] for a in alignmentsOfOneQuery)

def clumpSortKey(alignmentsPerQuery, clump):
    k = min(querySortKey(alignmentsPerQuery[i]) for i, isFlipped in clump)
    return -len(clump), k

def isInsertInRevStrand(insertSeqName, alignmentsOfOneQuery):
    for j, y in enumerate(alignmentsOfOneQuery):
        if j:
            x = alignmentsOfOneQuery[j - 1]
            if x[3] == insertSeqName and y[3] != insertSeqName:
                return y[8]
            if x[3] != insertSeqName and y[3] == insertSeqName:
                return x[8]

def isFlipQryStrand(opts, alignmentsOfOneQuery):
    if opts.insert:
        return isInsertInRevStrand(opts.insert, alignmentsOfOneQuery)
    return alignmentsOfOneQuery[0][8] and alignmentsOfOneQuery[-1][8]

def insertSortKey(opts, alignmentsPerQuery, clump):
    shift = 50
    seqName = opts.insert
    qryNum, isFlipped = clump[0]
    alns = alignmentsPerQuery[qryNum]
    for j, y in enumerate(alns):
        if j:
            x = alns[j - 1]
            xSeq, xBeg, xEnd = x[3:6]
            ySeq, yBeg, yEnd = y[3:6]
            if xSeq == seqName and ySeq != seqName:
                pos = yEnd - shift if y[8] else yBeg + shift
                return ySeq, pos
            if xSeq != seqName and ySeq == seqName:
                pos = xBeg + shift if x[8] else xEnd - shift
                return xSeq, pos

def isMergedGroupName(qryName):
    return re.match(r"(group|merged?)\d+-", qryName)

def groupIdFromName(mergedGroupName):
    return re.search(r"\d+", mergedGroupName).group()

def groupNameSortKey(mergedGroupName):
    return int(groupIdFromName(mergedGroupName))

def groupSortKey(queryNames, clump):
    return min(groupNameSortKey(queryNames[qryNum]) for qryNum, junk in clump)

def alignmentsInGenomeOrder(alignmentsPerQuery):
    logging.info("sorting...")
    alignments = chain.from_iterable(alignmentsPerQuery)
    return sorted(alignments, key=itemgetter(3, 4))

def jumpsPerChromosomeStrands(alignmentsPerQuery):
    jumpsDict = collections.defaultdict(list)
    for alignmentsOfOneQuery in alignmentsPerQuery:
        x = None
        for y in alignmentsOfOneQuery:
            yRefName = y[3]
            yStrand = y[8]
            if x:
                emptyListOfOverlaps = y[7]
                fwdKey = xStrand, xRefName, yStrand, yRefName
                revKey = not yStrand, yRefName, not xStrand, xRefName
                if fwdKey <= revKey:
                    jumpsDict[fwdKey].append((x, y, emptyListOfOverlaps))
                if revKey <= fwdKey:
                    jumpsDict[revKey].append((y, x, emptyListOfOverlaps))
            x, xRefName, xStrand = y, yRefName, yStrand
    return jumpsDict

def jumpSortKey(jump):
    return jump[0][4]

def alignmentsOfCoveredQueries(opts, alignmentsPerQuery):
    jumpsDict = jumpsPerChromosomeStrands(alignmentsPerQuery)
    logging.info("finding shared jumps...")
    for jumpsList in jumpsDict.values():
        jumpsList.sort(key=jumpSortKey)
        addSharedJumps(opts, jumpsList)
    while True:
        logging.info("excluding...")
        qryNums = set(i[0][0] for i in alignmentsPerQuery)
        alignmentsPerQuery = [i for i in alignmentsPerQuery
                              if isAllJumpsSupported(opts, qryNums, i)]
        if len(alignmentsPerQuery) == len(qryNums):
            break
        logging.info("queries: " + str(len(alignmentsPerQuery)))
    delOverlaps(chain.from_iterable(alignmentsPerQuery))
    return alignmentsPerQuery

def newStrand(strand, isFlipped):
    return "-+"[isFlipped == (strand == "-")]

def qryNameWithStrand(qryName, isFlipped):
    qryNameEnd = qryName[-1]
    isAddChar = qryNameEnd not in "+-"
    qryBase = qryName if isAddChar else qryName[:-1]
    newQryName = qryBase + newStrand(qryNameEnd, isFlipped)
    return newQryName, isAddChar

def printMaf(lines, isFlipped):
    lines = [re.split(r"(\s+)", i, 5) for i in lines]
    sLines = [i for i in lines if i[0] == "s"]
    qryLine = sLines[-1]
    qryName = qryLine[2]
    newQryName, isAddChar = qryNameWithStrand(qryName, isFlipped)
    qryLine[8] = newStrand(qryLine[8], isFlipped)
    sLineCount = 0
    for line in lines:
        if line[0] in "sq":
            if line[0] == "s":
                sLineCount += 1
            if sLineCount == len(sLines) and line[2] == qryName:
                line[2] = newQryName
            elif isAddChar:
                line[2] += " "
        elif line[0] == "p":
            if isAddChar:
                line[0] += " "
        print("".join(line), end="")
    print()

def printTab(fields, isFlipped):
    qryName = fields[6]
    newQryName, isAddChar = qryNameWithStrand(qryName, isFlipped)
    fields[6] = newQryName
    fields[9] = newStrand(fields[9], isFlipped)
    print(*fields, sep="\t")

def printShrunk(alignmentsOfOneQuery):
    oldRefName = None
    qryInc = 0
    for i in alignmentsOfOneQuery:
        qryNum, qryBeg, qryEnd, refName, refBeg, refEnd = i[:6]
        if qryBeg < 0:
            qryBeg, qryEnd = -qryEnd, -qryBeg
            refBeg, refEnd = -refEnd, -refBeg
        qryLen = qryEnd - qryBeg
        refLen = refEnd - refBeg
        refLenInc = refLen - qryLen
        if oldRefName:
            qryInc = qryBeg - oldQryEnd
        if refName != oldRefName:
            print(qryInc, qryLen, refBeg, refLenInc, refName, sep="\t")
        else:
            refInc = refBeg - oldRefEnd
            print(qryInc, qryLen, refInc, refLenInc, sep="\t")
        oldRefName = refName
        oldQryEnd = qryEnd
        oldRefEnd = refEnd
    print()

def printAlignments(opts, alnsPerKeptQuery, alignmentTextsPerQuery):
    for i in alnsPerKeptQuery:
        if opts.shrink:
            printShrunk(i)
        else:
            qryNum = i[0][0]
            for t in alignmentTextsPerQuery[qryNum]:
                separator = "" if t[0][0] == "a" else "\t"
                print(*t, sep=separator)

def alignmentsPerKeptQuery(opts, dataPerQry, alignments, ngFileNames):
    logging.info("queries: " + str(len(dataPerQry)))
    numOfNgSeqs = 4096  # start small: try to avoid too many overlaps
    for f in ngFileNames:
        ngIn = [i[0] for i in alignmentsPerRearrangedQuerySequence(opts, [f])]
        ngAlnsPerQry = iter(ngIn)
        ngAlignments = list(chain.from_iterable(ngIn))
        while True:
            ngAlns = alignmentsInGenomeOrder(islice(ngAlnsPerQry, numOfNgSeqs))
            if not ngAlns: break
            okAlns = alignmentsInGenomeOrder(i[0] for i in dataPerQry)
            addNgOverlaps(okAlns, ngAlns)
            logging.info("excluding...")
            dataPerQry = [i for i in dataPerQry
                          if isNoSharedRearrangement(opts, alignments,
                                                     ngAlignments, i[0], i[4])]
            logging.info("queries: " + str(len(dataPerQry)))
            delOverlaps(okAlns)
            numOfNgSeqs = min(numOfNgSeqs * 16, sys.maxsize)
    return [i[0] for i in dataPerQry]

def clumpsOfClumps(alignments, alignmentsPerQuery, clumps):
    links = set(linksBetweenClumps(alignments, alignmentsPerQuery, clumps))
    adjacencyList = adjacencyListFromLinks(len(clumps), links)
    isRevStrand = [False] * len(clumps)
    clumpPriorities = range(len(clumps))
    for i in connectedComponents(adjacencyList, clumpPriorities, isRevStrand):
        yield [clumps[clumpNum] for clumpNum, isFlipped in i]

def namedClumps(queryNames, isEachQueryOneMergedGroup, clumps):
    if isEachQueryOneMergedGroup:
        for clump in clumps:
            name = "merge" + "_".join(groupIdFromName(queryNames[qryNum])
                                      for qryNum, junk in clump)
            yield name, clump
    else:
        for i, clump in enumerate(clumps):
            name = "group{0}-{1}".format(i + 1, len(clump))
            yield name, clump

def wantedClumps(minNumOfFiles, fileNumPerQuery, clumps):
    for i in clumps:
        clumpName, clump = i
        fileNums = set(fileNumPerQuery[qryNum] for qryNum, isFlipped in clump)
        if len(fileNums) >= minNumOfFiles:
            yield i

def flippedAlignment(alignment, isFlipped):
    qryNum, qryBeg, qryEnd, refName, refBeg, refEnd = alignment[:6]
    # use reverse strand of query if isFlipped, else forward strand:
    if (qryBeg < 0) != isFlipped:
        qryBeg, qryEnd = -qryEnd, -qryBeg
        refBeg, refEnd = -refEnd, -refBeg
    return qryNum, qryBeg, qryEnd, refName, refBeg, refEnd

def rangeText(refRange):
    refName, refBeg, refEnd = refRange
    sign = ">" if refBeg < 0 else "<"
    return "{0}:{1}{2}{3}".format(refName, abs(refBeg), sign, abs(refEnd))

def isJoinableAlignments(opts, x, y):
    if refNameAndStrand(x) != refNameAndStrand(y):
        return False
    if y[4] >= x[5] + opts.min_gap:
        return False
    if y[1] >= x[2] + opts.min_gap:
        return False
    if x[5] >= y[4] + opts.min_rev or x[5] >= y[5]:
        return False
    return True

def refRangesFromFlippedAlns(opts, flippedAlns):
    for j, y in enumerate(flippedAlns):
        if j:
            x = flippedAlns[j - 1]
            if isJoinableAlignments(opts, x, y):
                refEnd = y[5]
                continue
            yield refName, refBeg, refEnd
        refName, refBeg, refEnd = y[3:6]
    yield refName, refBeg, refEnd

def qrySummary(opts, queryNames, alignmentsPerQuery, qryNum, isFlipped):
    qryName = queryNames[qryNum]
    newQryName, isAddChar = qryNameWithStrand(qryName, isFlipped)
    alns = alignmentsPerQuery[qryNum]
    flippedAlns = sorted(flippedAlignment(i, isFlipped) for i in alns)
    refRanges = refRangesFromFlippedAlns(opts, flippedAlns)
    texts = [rangeText(i) for i in refRanges]
    return newQryName, texts

def main(opts, args):
    logLevel = logging.INFO if opts.verbose else logging.WARNING
    logging.basicConfig(format="%(filename)s: %(message)s", level=logLevel)

    if not args:
        args.append("-")
    if ":" not in args:
        args.append(":")
    colonPos = args.index(":")
    colonEnd = colonPos + 1

    okInput = list(alignmentsPerRearrangedQuerySequence(opts, args[:colonPos]))
    alignmentsPerQuery = [i[0] for i in okInput]
    alignmentTextsPerQuery = [i[1] for i in okInput]
    fileNumPerQuery = [i[2] for i in okInput]
    queryNames = [i[3] for i in okInput]
    alignments = list(chain.from_iterable(alignmentsPerQuery))

    alnsPerKeptQuery = alignmentsPerKeptQuery(opts, okInput, alignments,
                                              args[colonEnd:])

    if opts.min_cov:
        alnsPerKeptQuery = alignmentsOfCoveredQueries(opts, alnsPerKeptQuery)

    print("#", os.path.basename(sys.argv[0]), *sys.argv[1:])
    print()

    if opts.min_seqs < 1:
        return printAlignments(opts, alnsPerKeptQuery, alignmentTextsPerQuery)

    addOverlaps(alignmentsInGenomeOrder(alnsPerKeptQuery))
    links = linksBetweenQueries(opts, alignments, alnsPerKeptQuery)
    adjacencyList = adjacencyListFromLinks(len(alignmentsPerQuery), links)

    logging.info("grouping...")
    isRevStrand = [isFlipQryStrand(opts, i) for i in alignmentsPerQuery]
    qryPriorities = [(-len(i), -alignedQueryLength(j))
                     for i, j in zip(adjacencyList, alignmentsPerQuery)]
    allClumps = connectedComponents(adjacencyList, qryPriorities, isRevStrand)
    clumps = [i for i in allClumps if len(i) >= opts.min_seqs]
    keptQueryNums = set(i[0][0] for i in alnsPerKeptQuery)
    clumps = [i for i in clumps if i[0][0] in keptQueryNums]
    isEachQueryOneMergedGroup = all(map(isMergedGroupName, queryNames))
    sortKey = functools.partial(clumpSortKey, alignmentsPerQuery)
    if opts.insert:
        sortKey = functools.partial(insertSortKey, opts, alignmentsPerQuery)
    if isEachQueryOneMergedGroup:
        sortKey = functools.partial(groupSortKey, queryNames)
    clumps.sort(key=sortKey)
    clumpClumps = clumpsOfClumps(alignments, alignmentsPerQuery, clumps)
    clumps = chain.from_iterable(clumpClumps)
    clumps = namedClumps(queryNames, isEachQueryOneMergedGroup, clumps)
    goodClumps = list(wantedClumps(colonPos, fileNumPerQuery, clumps))

    if not opts.shrink:
        for clumpName, clump in goodClumps:
            print("#", clumpName)
            s = [qrySummary(opts, queryNames, alignmentsPerQuery, *i)
                 for i in clump]
            width = max(len(newQryName) for newQryName, ranges in s)
            for newQryName, ranges in s:
                paddedName = format(newQryName, str(width))
                para = " ".join([paddedName] + ranges)
                wrapWidth = opts.width if opts.width else len(para) + 2
                print(textwrap.fill(para, wrapWidth, break_long_words=False,
                                    break_on_hyphens=False,
                                    initial_indent="# ",
                                    subsequent_indent="#  "))
            print()

    for clumpName, clump in goodClumps:
        if opts.shrink:
            for qryNum, isFlipped in clump:
                printShrunk(alignmentsPerQuery[qryNum])
            print()
        else:
            print("# PART", clumpName)
            print()
            for qryNum, isFlipped in clump:
                for t in alignmentTextsPerQuery[qryNum]:
                    if t[0][0] == "a":
                        printMaf(t, isFlipped)
                    elif len(t) > 11:
                        printTab(t, isFlipped)
                    else:  # shrunk format:
                        print(*t, sep="\t")
                print()

    sys.stdout.flush()
    sys.stderr.write("number of groups: {0}\n".format(len(goodClumps)))

if __name__ == "__main__":
    signal.signal(signal.SIGPIPE, signal.SIG_DFL)  # avoid silly error message
    usage = "%prog [options] case-file(s) [: control-file(s)]"
    descr = "Find rearranged query sequences in query-to-reference alignments."
    op = optparse.OptionParser(usage=usage, description=descr)
    op.add_option("-s", "--min-seqs", type="int", default=2, metavar="N",
                  help="minimum query sequences per group (default=%default)")
    op.add_option("-c", "--min-cov", type="int", metavar="N", help=
                  "omit any query with any rearrangement shared by < N "
                  "other queries (default: 1 if s>1, else 0)")
    op.add_option("-t", "--types", metavar="LETTERS", default="CSNG", help=
                  "rearrangement types: C=inter-chromosome, S=inter-strand, "
                  "N=non-colinear, G=big gap (default=%default)")
    op.add_option("-g", "--min-gap", type="float", default=10000, metavar="BP",
                  help='minimum forward jump in the reference sequence '
                  'counted as a "big gap" (default=%default)')
    op.add_option("-r", "--min-rev", type="float", default=1000, metavar="BP",
                  help='minimum reverse jump in the reference sequence '
                  'counted as "non-colinear" (default=%default)')
    op.add_option("-f", "--filter", type="int", default=1, metavar="N",
                  help='discard case reads sharing any (0) or "strongest" '
                  '(1) rearrangements with control reads (default=%default)')
    op.add_option("-d", "--max-diff", type="float", default=500, metavar="BP",
                  help="maximum query-length difference for "
                  "shared rearrangement (default=%default)")
    op.add_option("-m", "--max-mismap", type="float", default=1.0,
                  metavar="PROB", help="discard any alignment with "
                  "mismap probability > PROB (default=%default)")
    op.add_option("--insert", metavar="NAME",
                  help="find insertions of the sequence with this name")
    op.add_option("--shrink", action="count", help="shrink the output")
    op.add_option("-v", "--verbose", action="count",
                  help="show progress messages")
    op.add_option("-w", "--width", type="int", default=79, metavar="W", help=
                  "line-wrap width of group summary lines (default=%default)")
    opts, args = op.parse_args()
    if opts.min_cov is None:
        opts.min_cov = 1 if opts.min_seqs > 1 else 0
    main(opts, args)