def findRegions(tranList):
    # Find breakpoints where coverage by exons changes.

    # Why are we doing this? See the note in the Transcript class
    # definition below.

    breaks = list()
    for tranIx, tran in enumerate(tranList):
        for exon in tran.exons:
            breaks.append([exon.start, 0, tranIx, tran.name, exon.name])
            breaks.append([exon.end, 1, tranIx, tran.name, exon.name])

    breaks.sort(key=lambda x: x[0])
    curPos = breaks[0][0]
    curTranSet = set()
    region = 0

    for ix in xrange(len(breaks)):
        posit, flag, tranIx, tranName, exonName = breaks[ix]
        if posit > curPos + MIN_REGION_SIZE:             # this is a new region
            if len(curTranSet) > 0:
                for ix in curTranSet:
                    tranList[ix].regions.add(region)     # update set of regions hit by this transcript
                region += 1
            curPos = posit
        if flag == 0:                                    # exon start
            curTranSet.add(tranIx)
        else:                                            # exon end
            curTranSet.remove(tranIx)
    logger.debug('found %d regions' % region)
    return

    def __init__(self, filename, CCSDir=None):

        logger.debug("creating BasFile object")

        self._filename = filename
        self._CCSDir = CCSDir
        self._infile = h5py.File(filename, "r")
        self._top = self._infile
        self._baxfile = list()
        self._coords = None

        self._consensusIndex = None
        self._consPassIndex = None

        if "MultiPart" not in self._top:  # if this is an old-style bas file
            bf = BaxFile(filename)  # file will contain its own CCSdata
            self._baxfile.append(bf)  # only one file (this one) in the list

        else:  # else it's an index to a set of bax files

            h5Dir = os.path.dirname(os.path.abspath(self._filename))

            for baxfileName in self._top["MultiPart/Parts"]:  # for each bax file

                fqBaxfileName = os.path.join(h5Dir, baxfileName)  # fq = fully qualified
                bf = BaxFile(fqBaxfileName, CCSDir=CCSDir)
                self._baxfile.append(bf)  # add file to list

        self.fillCombinedFields()  # need to compute this first, we'll need it later
        self.fillZMWIndexes()
        self.fillMovieName()
        self.fillRegionIndexes()

    def fillConsensusIndexes(self):
        """Compute _consPassIndex and _consensusIndex. Called only when those arrays are needed."""

        self._consPassIndex = [None] * (self._maxZMW + 1)
        self._consensusIndex = [None] * (self._maxZMW + 1)

        for bf in self._baxfile:

            bf.ZMWSanityClause()  # sanity check the consensus datasets

            ix = 0
            passIndex = 0
            consIndex = 0
            numPasses = bf._consPasses["NumPasses"]
            numEvent = bf._consZMW["NumEvent"]

            for hole in bf.holeNumbers():
                self._consPassIndex[hole] = passIndex
                passIndex += numPasses[ix]
                self._consensusIndex[hole] = consIndex
                consIndex += numEvent[ix]
                ix += 1

            logger.debug("%s processed %d consensus passes" % (bf.shortName(), passIndex))
            logger.debug("%s processed %d consensus basecalls" % (bf.shortName(), consIndex))

    def basecallIndex (self):
        '''Create and cache a table by hole number of starting indexes into PulseData/BaseCalls.'''

        if self._basecallIndex == None:

            logger.debug("creating basecall index")
            
            numEvent   = self._ZMW["NumEvent"]
            holeNumber = self._ZMW["HoleNumber"]
            numZ       = self.numZMWs()
            index      = 0
            self._basecallIndex = [0] * numZ

            # The loop below includes a check that the HoleNumbers run
            # from 0 to N-1. I.e., HoleNumber[ix] == ix. Otherwise,
            # there is no way to get from the hole number of a region
            # back to the ZMW entry (other than searching for it).

            for ix in xrange(numZ):
                
                if holeNumber[ix] != ix:
                    raise RuntimeError("Hole number != index at %d" % ix)

                self._basecallIndex[ix] = index
                index += numEvent[ix]

            logger.debug("processed %d basecalls" % index)

        return self._basecallIndex

    def cellCoords (self):
        '''Find and cache the minimum and maximum X/Y coordinates on the SMRTcell'''

        if self._coords is None:

            logger.debug("finding SMRTcell coordinates")

            numZ       = self.numZMWs()
            holeXY = self._ZMW["HoleXY"]

            minX = maxX = holeXY[0,0]
            minY = maxY = holeXY[0,1]

            for ix in xrange(numZ):

                x,y = holeXY[ix,:]

                if x < minX:
                    minX = x
                elif x > maxX:
                    maxX = x

                if y < minY:
                    minY = y
                elif y > maxY:
                    maxY = y

            self._coords = (minX, maxX, minY, maxY)

            logger.debug("SMRTcell is (%d,%d,%d,%d)" % (minX, maxX, minY, maxY))

        return self._coords

    def cellCoords(self):
        """Find and cache the minimum and maximum X/Y coordinates on the SMRTcell"""

        if self._coords is None:

            logger.debug("finding SMRTcell coordinates")

            minX = 0
            maxX = 0
            minY = 0
            maxY = 0

            for bf in self._baxfile:

                holeXY = bf._ZMW["HoleXY"]

                minX = min(minX, min(holeXY[:, 0]))
                maxX = max(maxX, max(holeXY[:, 0]))
                minY = min(minY, min(holeXY[:, 1]))
                maxY = max(maxY, max(holeXY[:, 1]))

            self._coords = (minX, maxX, minY, maxY)

            logger.debug("SMRTcell is (%d,%d,%d,%d)" % (minX, maxX, minY, maxY))

        return self._coords

    def findCCSFile(self):
        """Given a directory to look in, find the ccs.h5 file that contains consensus reads for this bax file."""

        self._hasConsensus = False  # until proven otherwise

        if "PulseData/ConsensusBaseCalls" in self._top:  # if this is an older bax file, in contains its own CCS data

            self._consBasecalls = self._top["PulseData/ConsensusBaseCalls"]
            self._consZMW = self._top["PulseData/ConsensusBaseCalls/ZMW"]
            self._consPasses = self._top["PulseData/ConsensusBaseCalls/Passes"]
            self._hasConsensus = True

        elif self._CCSDir is not None:

            CCSFilename = os.path.basename(self._filename).replace("bax", "ccs")
            fqCCSFilename = os.path.join(self._CCSDir, CCSFilename)

            if os.path.exists(fqCCSFilename):

                self._CCSFile = h5py.File(fqCCSFilename, "r")
                self._consBasecalls = self._CCSFile["PulseData/ConsensusBaseCalls"]
                self._consZMW = self._CCSFile["PulseData/ConsensusBaseCalls/ZMW"]
                self._consPasses = self._CCSFile["PulseData/ConsensusBaseCalls/Passes"]
                self._hasConsensus = True
                logger.debug("BaxFile %s found CCS file %s" % (self._shortName, fqCCSFilename))

            else:
                logger.warning("%s: no CCS file found corresponding to %s" % (self._shortName, self._filename))

        else:
            logger.info("BaxFile %s does not contain CCS data (rel 2.1.0 and later). Use --ccs" % self._shortName)

    def _fillRegionTables (self):
        '''Create and cache tables by hole number of starting and HQ indexes into PulseData/Regions.'''

        if self._regionIndex == None:

            logger.debug("creating region index")
            
            self._regionIndex = [0] * self.numZMWs()
            self._HQIndex     = [0] * self.numZMWs()

            regions  = self._regions
            index    = 0
            lastHole = -1        # init to non-matching value

            for line in regions:

                hole, regionType = line[0:2]

                if hole != lastHole:                    # start of new hole?
                    self._regionIndex[hole] = index
                    lastHole = hole

                if regionType == 2:                     # HQ region for this hole?
                    self._HQIndex[hole] = index

                index += 1
                    
            logger.debug("processed %d regions" % index)

        return self._regionIndex

    def __init__ (self, name):

        self.name = name
        self.handle = open (name, 'w')
        self.lastPos = dict()

        logger.debug('opened %s' % name)

    def toPickle (self, filename):

        pickHandle = open (filename, 'w')
        pk = pickle.Pickler (pickHandle, pickle.HIGHEST_PROTOCOL)
        pk.dump (self)
        pickHandle.close()

        logger.debug('wrote annotation data to pickle file %s' % filename)

        return

    def __init__ (self, fileName):

        logger.debug("creating CmpFile object for %s" % (fileName))

        self._fileName = fileName
        self._infile   = h5py.File (fileName, 'r')
####        self._top      = h5py.Group (self._infile, '/')
        self._top      = self._infile         # h5py 2.0.1 change!

        return

    def fromPickle (filename):
        '''Create an AnnotationList object from a pickle file (alternative to __init__).'''

        logger.debug('reading annotations in pickle format from %s' % filename)

        handle = open (filename, 'r')
        pk = pickle.Unpickler (handle)
        annotList = pk.load()
        handle.close()

        return annotList

    def fromPickle (filename):
        '''Create a Reference object from a pickle file (alternative to __init__).'''

        logger.debug('reading reference in pickle format from %s' % filename)

        handle = open (filename, 'r')
        pk = pickle.Unpickler (handle)
        ref = pk.load()
        handle.close()

        return ref

    def fromPickle (filename):
        '''Create a ClusterDict object from a pickle file (alternative to __init__).'''

        handle = open (filename, 'r')
        pk = pickle.Unpickler (handle)
        clusterDict = pk.load()
        handle.close()

        logger.debug('read %d clusters in pickle format from %s' % (len(clusterDict), filename))

        return clusterDict

def submitFinalJobs (opt, chunkList):

    chunkFiles = ['%s \\\n' % chk.trimmedChunkName for chk in chunkList]

    sh = list()
    sh.append('#!/bin/bash\n\n')
    sh.append('set -o errexit\n')
    sh.append('set -o nounset\n\n')

    sh.append('cat \\\n')
    sh.extend(chunkFiles)
    sh.append(' > %s\n' % opt.output)

    if opt.report is not None:
        reportFiles = ['%s \\\n' % chk.reportChunkName for chk in chunkList]
        sh.append('\ncat \\\n')
        sh.extend(reportFiles)
        sh.append(' > %s\n' % opt.report)

    finalScriptName = '%s/trim_final.sh' % opt.tmpdir
    handle =  open (finalScriptName, 'w')
    handle.writelines (sh)
    handle.close()

    deps = ':'.join ([chk.jobno for chk in chunkList])

    cmd = list()
    cmd.append('qsub')
    cmd.append('-N trim_final')       # job name
    cmd.append('-o trim_final.out')   # output file
    cmd.append('-j oe')               # combine stdout and stderr
    cmd.append('-l nodes=1:ppn=1,walltime=4:00:00')    # resources required
    cmd.append('-d . ')               # working directory (strangely, ./ is not the default)
    cmd.append('-r n')                # do NOT attempt to restart on failure
    cmd.append('-V')                  # export all environment variables to job
    cmd.append('-W umask=0002')       # make logs rw-rw-r--
    cmd.append('-m n')                # don't send any mail
    cmd.append('-W depend=afterok:%s' % deps)
    cmd.append(finalScriptName)       # script to run

    command = ' '.join(cmd)
    logger.debug ('running %s' % command)
    
    popen_file = os.popen(command)
    response = popen_file.read().strip()
    rc = popen_file.close()
    if rc is not None:
        logger.error('command failed, rc=%d' % rc)
        raise RuntimeError

    logger.debug ('jobno is %s' % response)

    return response

def plotPolyAs (tranList, blocks):
    '''Add start/stop codons to plot.'''

    for tran in tranList:
        if tran.annot:                                        # only annotations know about polyAs
            for exon in tran.exons:
                if hasattr (exon, 'polyAs'):
                    for start, end, howmany in exon.polyAs:
                        plotA (tran, start, blocks)
                        logger.debug ('%s: %9d' % (exon.name, start))

    return

    def toPickle (self, filename):

        self.geneDict =  None         # no need to pickle this, it can be recreated

        pickHandle = open (filename, 'w')
        pk = pickle.Pickler (pickHandle, pickle.HIGHEST_PROTOCOL)
        pk.dump (self)
        pickHandle.close()

        logger.debug('wrote %d clusters to pickle file %s' % (len(self.clusterDict), filename))

        return

    def makeRef (self):
        '''Invoke bowtie-build on the fasta file.'''

        if not self.handle.closed:
            self.close()

        command = '%s %s %s > %s.out 2>&1' % (BOWTIE_BUILD, self.name, self.name, self.name)
        logger.debug(command)

        buildOut = os.popen (command)      # this should return nothing, since we've redirected the output
        rc = buildOut.close()
        if rc is not None:
            raise RuntimeError ('bowtie2-build failed: %d' % rc)

def main ():

    logger.debug('version %s starting' % VERSION)

    opt, args = getParms()

    makeTempDir (opt.tmpdir)

    nSeqs = countSeqs (opt.input)
    logger.debug('%s contains %d sequences' % (opt.input, nSeqs))

    seqsPerJob = (nSeqs + opt.njobs - 1) / opt.njobs
    logger.debug('each of %d jobs will process %d sequences' % (opt.njobs, seqsPerJob))

    chunkList = makeFastaChunks (opt, nSeqs, seqsPerJob)

    for chunk in chunkList:
        chunk.makeScript()
        chunk.submitScript()

    submitFinalJobs (opt, chunkList)

    logger.debug('finished')

    return

    def getGeneDict (self):
        '''Create and cache dict: key=gene name, value=Annotation object for gene.'''

        if self.geneDict is None:

            logger.debug('creating gene name lookup table')

            self.geneDict = dict()

            for chr in self.chromosomes():
                for gene in self.annot[chr].getChildren():
                    self.geneDict.setdefault(gene.name, []).append(gene)

        return self.geneDict