class TimeKeeper:
    def __init__(self):
        self.logger = logging.getLogger('timestamp')
        self.startTime = time.time()
        self.lastLogTime = self.startTime

    def startTimer(self):
        """Restart the timer"""
        self.startTime = time.time()
        self.lastLogTime = self.startTime

    def getTimeStamp(self):
        """Make a time stamp"""
        now = time.time()
        ret_str = "{ Current stage: %s || Total: %s }" % (self.secondsToStr(now - self.lastLogTime), self.secondsToStr(now - self.startTime))
        self.lastLogTime = now
        return ret_str

    def printTimeStamp(self):
        self.logger.info(self.getTimeStamp())

    def secondsToStr(self, t):
        rediv = lambda ll, b: list(divmod(ll[0], b)) + ll[1:]
        return "%d:%02d:%02d.%03d" % tuple(reduce(rediv, [[t * 1000, ], 1000, 60, 60])) class MarkerSetParser():
    """Parse marker set file."""

    def __init__(self, threads=1):
        self.logger = logging.getLogger('timestamp')
        self.numThreads = threads

    def getMarkerSets(self, outDir, binIds, markerFile, excludeMarkersFile=None):
        """Determine marker set for each bin."""

        # determine type of marker set file
        markerFileType = self.markerFileType(markerFile)

        # get marker set for each bin
        binIdToBinMarkerSets = {}

        if markerFileType == BinMarkerSets.TAXONOMIC_MARKER_SET:
            binMarkerSets = self.parseTaxonomicMarkerSetFile(markerFile)

            for binId in binIds:
                binIdToBinMarkerSets[binId] = binMarkerSets
        elif markerFileType == BinMarkerSets.TREE_MARKER_SET:
            binIdToBinMarkerSets = self.parseLineageMarkerSetFile(markerFile)
        else:
            markers = [set()]
            modelParser = HmmModelParser(markerFile)
            for model in modelParser.parse():
                markers[0].add(model.acc)
            markerSet = MarkerSet(0, "N/A", -1, markers)

            for binId in binIds:
                binMarkerSets = BinMarkerSets(binId, BinMarkerSets.HMM_MODELS_SET)
                binMarkerSets.addMarkerSet(markerSet)
                binIdToBinMarkerSets[binId] = binMarkerSets

        # remove marker genes specified by user or marker for exclusion
        markersToExclude = set()
        if excludeMarkersFile:
            markersToExclude = self.readExcludeMarkersFile(excludeMarkersFile)

        markersToExclude.update(DefaultValues.MARKERS_TO_EXCLUDE)
        for binId, binMarkerSet in binIdToBinMarkerSets.items():
            binMarkerSet.removeMarkers(markersToExclude)

        return binIdToBinMarkerSets

    def readExcludeMarkersFile(self, excludeMarkersFile):
        """Parse file specifying markers to exclude."""
        markersToExclude = set()
        for line in open(excludeMarkersFile):
            if line[0] == '#':
                continue

            marker = line.strip()
            markersToExclude.add(marker)

        return markersToExclude

    def createHmmModels(self, outDir, binIds, markerFile):
        """Create HMM model for each bins marker set."""

        # determine type of marker set file
        markerFileType = self.markerFileType(markerFile)

        # get HMM file for each bin
        binIdToModels = {}
        if markerFileType == BinMarkerSets.TAXONOMIC_MARKER_SET:
            hmmModelFile = self.createHmmModelFile(list(binIds.keys())[0], markerFile)

            modelParser = HmmModelParser(hmmModelFile)
            models = modelParser.models()
            for binId in binIds:
                binIdToModels[binId] = models

            os.remove(hmmModelFile)
        elif markerFileType == BinMarkerSets.TREE_MARKER_SET:
            binIdToModels = self.__createLineageHmmModels(binIds, markerFile)
        else:
            modelParser = HmmModelParser(markerFile)
            models = modelParser.models()
            for binId in binIds:
                binIdToModels[binId] = models

        return binIdToModels

    def createHmmModelFile(self, binId, markerFile):
        """Create HMM file for from a bin's marker set."""

        # determine type of marker set file
        markerFileType = self.markerFileType(markerFile)

        # create HMM file
        hmmModelFile = os.path.join(tempfile.gettempdir(), str(uuid.uuid4()))
        if markerFileType == BinMarkerSets.TAXONOMIC_MARKER_SET:
            binMarkerSets = self.parseTaxonomicMarkerSetFile(markerFile)
            self.__createMarkerHMMs(binMarkerSets, hmmModelFile, bReportProgress=False)
        elif markerFileType == BinMarkerSets.TREE_MARKER_SET:
            binIdToBinMarkerSets = self.parseLineageMarkerSetFile(markerFile)
            self.__createMarkerHMMs(binIdToBinMarkerSets[binId], hmmModelFile, bReportProgress=False)
        else:
            shutil.copyfile(markerFile, hmmModelFile)

        return hmmModelFile

    def __createLineageHmmModels(self, binIds, markerFile):
        """Create lineage-specific HMMs for each bin."""

        self.logger.info('Extracting lineage-specific HMMs with %d threads:' % self.numThreads)

        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        for binId in binIds:
            workerQueue.put(binId)

        for _ in range(self.numThreads):
            workerQueue.put(None)

        binIdToModels = mp.Manager().dict()

        try:
            calcProc = [mp.Process(target=self.__fetchModelInfo, args=(binIdToModels, markerFile, workerQueue, writerQueue)) for _ in range(self.numThreads)]
            writeProc = mp.Process(target=self.__reportFetchProgress, args=(len(binIds), writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put(None)
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

        # create a standard dictionary from the managed dictionary
        d = {}
        for binId in binIdToModels.keys():
            d[binId] = binIdToModels[binId]

        return d

    def __fetchModelInfo(self, binIdToModels, markerFile, queueIn, queueOut):
        """Fetch HMM."""
        while True:
            binId = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            hmmModelFile = self.createHmmModelFile(binId, markerFile)

            modelParser = HmmModelParser(hmmModelFile)
            binIdToModels[binId] = modelParser.models()

            os.remove(hmmModelFile)

            queueOut.put(binId)

    def __reportFetchProgress(self, numBins, queueIn):
        """Report progress of extracted HMMs."""

        numProcessedBins = 0
        if self.logger.getEffectiveLevel() <= logging.INFO:
                statusStr = '    Finished extracting HMMs for %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        while True:
            binId = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            if self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedBins += 1
                statusStr = '    Finished extracting HMMs for %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

    def markerFileType(self, markerFile):
        """Determine type of marker file."""
        with open(markerFile, 'r') as f:
            header = f.readline()

        if DefaultValues.TAXON_MARKER_FILE_HEADER in header:
            return BinMarkerSets.TAXONOMIC_MARKER_SET
        elif DefaultValues.LINEAGE_MARKER_FILE_HEADER in header:
            return BinMarkerSets.TREE_MARKER_SET
        elif 'HMMER3' in header:
            return BinMarkerSets.HMM_MODELS_SET
        else:
            self.logger.error('Unrecognized file type: ' + markerFile)
            sys.exit(1)

    def __createMarkerHMMs(self, binMarkerSet, outputFile, bReportProgress=True):
        """Create HMM file for markers."""

        # get list of marker genes
        markerGenes = binMarkerSet.getMarkerGenes()

        # get all genes from the same clan as any marker gene
        pfam = PFAM(DefaultValues.PFAM_CLAN_FILE)
        genesInSameClan = pfam.genesInSameClan(markerGenes)

        # extract marker genes along with all genes from the same clan
        allMarkers = markerGenes | genesInSameClan

        if bReportProgress:
            self.logger.info("There are %d genes in the marker set and %d genes from the same PFAM clan." % (len(markerGenes), len(genesInSameClan)))

        # create file with all model accession numbers
        keyFile = os.path.join(tempfile.gettempdir(), str(uuid.uuid4()))
        fout = open(keyFile, 'w')
        for modelAcc in allMarkers:
            fout.write(modelAcc + '\n')
        fout.close()

        # fetch specified models
        HF = HMMERRunner(mode='fetch')
        HF.fetch(DefaultValues.HMM_MODELS, keyFile, outputFile, bKeyFile=True)

        # index the HMM file
        if os.path.exists(outputFile + '.ssi'):
            os.remove(outputFile + '.ssi')
        HF.index(outputFile)

        # remove key file
        os.remove(keyFile)

    def parseTaxonomicMarkerSetFile(self, markerSetFile):
        """Parse marker set from a taxonomic-specific marker set file."""
        with open(markerSetFile) as f:
            f.readline()  # skip header

            binLine = f.readline()
            taxonId = binLine.split('\t')[0]
            binMarkerSets = BinMarkerSets(taxonId, BinMarkerSets.TAXONOMIC_MARKER_SET)
            binMarkerSets.read(binLine)

        return binMarkerSets

    def parseLineageMarkerSetFile(self, markerSetFile):
        """Parse marker sets from a lineage-specific marker set file."""

        # read all marker sets
        binIdToBinMarkerSets = {}
        with open(markerSetFile) as f:
            f.readline()  # skip header

            for line in f:
                lineSplit = line.split('\t')
                binId = lineSplit[0]

                binMarkerSets = BinMarkerSets(binId, BinMarkerSets.TREE_MARKER_SET)
                binMarkerSets.read(line)

                # determine selected marker set
                selectedMarkerSetMap = self.parseSelectedMarkerSetMap()
                binMarkerSets.setLineageSpecificSelectedMarkerSet(selectedMarkerSetMap)

                binIdToBinMarkerSets[binId] = binMarkerSets

        return binIdToBinMarkerSets

    def parseSelectedMarkerSetMap(self):
        selectedMarkerSetMap = {}
        for line in open(DefaultValues.SELECTED_MARKER_SETS):
            lineSplit = line.split('\t')
            internalID = lineSplit[0]
            selectedID = lineSplit[1].rstrip()

            selectedMarkerSetMap[internalID] = selectedID

        return selectedMarkerSetMap

    def writeBinModels(self, binIdToModels, filename):
        """Save HMM model info for each bin to file."""

        self.logger.info('Saving HMM info to file.')

        with gzip.open(filename, 'wb') as output:
            pickle.dump(binIdToModels, output, pickle.HIGHEST_PROTOCOL)

    def loadBinModels(self, filename):
        """Read HMM model info for each bin from file."""

        self.logger.info('Reading HMM info from file.')

        with gzip.open(filename, 'rb') as f:
            binIdToModels = pickle.load(f)

        return binIdToModels
 class BinMarkerSets():
    """A collection of one or more marker sets associated with a bin."""

    # type of marker set
    TAXONOMIC_MARKER_SET = 1
    TREE_MARKER_SET = 2
    HMM_MODELS_SET = 3

    def __init__(self, binId, markerSetType):
        self.logger = logging.getLogger('timestamp')
        self.markerSets = []
        self.binId = binId
        self.markerSetType = markerSetType
        self.selectedLinageSpecificMarkerSet = None

    def numMarkerSets(self):
        """Number of marker sets associated with bin."""
        return len(self.markerSets)

    def addMarkerSet(self, markerSet):
        """Add marker set to bin."""
        self.markerSets.append(markerSet)

    def markerSetIter(self):
        """Generator function for iterating over marker sets."""
        for markerSet in self.markerSets:
            yield markerSet

    def getMarkerGenes(self):
        """Get marker genes from all marker sets."""
        markerGenes = set()
        for ms in self.markerSets:
            markerGenes.update(ms.getMarkerGenes())

        return markerGenes

    def mostSpecificMarkerSet(self):
        return self.markerSets[0]

    def treeMarkerSet(self):
        pass

    def selectedMarkerSet(self):
        """Return the 'selected' marker set for this bin."""
        if self.markerSetType == self.TAXONOMIC_MARKER_SET:
            return self.mostSpecificMarkerSet()
        elif self.markerSetType == self.TREE_MARKER_SET:
            return self.selectedLinageSpecificMarkerSet
        else:
            # there should be a single marker set associate with this bin
            if len(self.markerSets) == 1:
                return self.markerSets[0]

            self.logger.error('Expect a single marker set to be associated with each bin.')
            sys.exit(1)

    def setLineageSpecificSelectedMarkerSet(self, selectedMarkerSetMap):
        uid = self.mostSpecificMarkerSet().UID

        selectedId = selectedMarkerSetMap[uid]
        self.selectedLinageSpecificMarkerSet = None
        
        while not self.selectedLinageSpecificMarkerSet:
            for ms in self.markerSets:
                if ms.UID == selectedId:
                    self.selectedLinageSpecificMarkerSet = ms
                    break
                
            if not self.selectedLinageSpecificMarkerSet:
                # This is a hack for the reduced tree. Since not all
                # marker sets are in the reduced tree it is possible the
                # selected marker set might not be avaliable. In this case,
                # we should move to the next suitable marker set. Ideally,
                # this could be avoided by just forcing in the selected
                # marker set.
                selectedId = selectedMarkerSetMap[selectedId]
            else:
                break

        if self.selectedLinageSpecificMarkerSet == None:
            # something has gone wrong
            self.logger.error('Failed to set a selected lineage-specific marker set.')
            sys.exit(1)

    def removeMarkers(self, markersToRemove):
        """Remove specified markers from all marker sets."""
        for markerSet in self.markerSets:
            curMarkerSet = markerSet.markerSet

            newMarkerSet = []
            for ms in curMarkerSet:
                newMS = ms - markersToRemove

                if len(newMS) != 0:
                    newMarkerSet.append(newMS)

            markerSet.markerSet = newMarkerSet

    def write(self, fout):
        """Write marker set to file."""
        fout.write(self.binId)
        fout.write('\t' + str(len(self.markerSets)))
        for ms in self.markerSets:
            fout.write('\t' + str(ms))
        fout.write('\n')

    def read(self, line):
        """Construct bin marker set data from line."""
        lineSplit = line.split('\t')
        numMarkerSets = int(lineSplit[1])
        for i in range(0, numMarkerSets):
            uid = lineSplit[i * 4 + 2]
            lineageStr = lineSplit[i * 4 + 3]
            numGenomes = int(lineSplit[i * 4 + 4])
            markerSet = eval(lineSplit[i * 4 + 5])
            self.markerSets.append(MarkerSet(uid, lineageStr, numGenomes, markerSet))
 class ResultsParser():
    """Parse output of Prodigal+HMMER run and derived statistics."""
    def __init__(self, binIdToModels):
        self.logger = logging.getLogger('timestamp')

        self.results = {}
        self.models = binIdToModels

    def analyseResults(self,
                       outDir,
                       binStatsFile,
                       hmmTableFile,
                       bIgnoreThresholds=False,
                       evalueThreshold=DefaultValues.E_VAL,
                       lengthThreshold=DefaultValues.LENGTH,
                       bSkipPseudoGeneCorrection=False,
                       bSkipAdjCorrection=False,
                       ):
        """Parse results in the output directory."""

        # read bin and sequence stats into dictionaries
        binStats = self.parseBinStats(outDir, binStatsFile)

        # get hits to each bin
        self.parseBinHits(outDir, hmmTableFile, bSkipAdjCorrection, bIgnoreThresholds, evalueThreshold, lengthThreshold, bSkipPseudoGeneCorrection, binStats)

        return binStats

    def cacheResults(self, outDir, binIdToBinMarkerSets, bIndividualMarkers):
        # cache critical results to file
        self.__writeBinStatsExt(outDir, binIdToBinMarkerSets, bIndividualMarkers)
        self.__writeMarkerGeneStats(outDir, binIdToBinMarkerSets, bIndividualMarkers)

    def parseBinHits(self, outDir,
                     hmmTableFile,
                     bSkipAdjCorrection=False,
                     bIgnoreThresholds=False,
                     evalueThreshold=DefaultValues.E_VAL,
                     lengthThreshold=DefaultValues.LENGTH,
                     bSkipPseudoGeneCorrection=False,
                     binStats=None):
        """ Parse HMM hits for each bin."""
        
        if not self.models:
            self.logger.error('Models must be parsed before identifying HMM hits.')
            sys.exit(1)

        self.logger.info('Parsing HMM hits to marker genes:')

        numBinsProcessed = 0
        for binId in self.models:
            if self.logger.getEffectiveLevel() <= logging.INFO:
                numBinsProcessed += 1
                statusStr = '    Finished parsing hits for %d of %d (%.2f%%) bins.' % (numBinsProcessed, len(self.models), float(numBinsProcessed) * 100 / len(self.models))
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

            if binStats != None:
                resultsManager = ResultsManager(binId, self.models[binId], bIgnoreThresholds, evalueThreshold, lengthThreshold, bSkipPseudoGeneCorrection, binStats[binId])
            elif binStats == None:
                resultsManager = ResultsManager(binId, self.models[binId], bIgnoreThresholds, evalueThreshold, lengthThreshold, bSkipPseudoGeneCorrection)
            else:
                self.logger.error('Invalid parameter settings for binStats and seqStats.')
                sys.exit(1)

            hmmerTableFile = os.path.join(outDir, 'bins', binId, hmmTableFile)
            self.parseHmmerResults(hmmerTableFile, resultsManager, bSkipAdjCorrection)
            self.results[binId] = resultsManager

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

    def __writeBinStatsExt(self, directory, binIdToBinMarkerSets, bIndividualMarkers):
        binStatsExtFile = os.path.join(directory, 'storage', DefaultValues.BIN_STATS_EXT_OUT)
        fout = open(binStatsExtFile, 'w')

        for binId in self.results:
            binStatsExt = self.results[binId].getSummary(binIdToBinMarkerSets[binId], bIndividualMarkers, outputFormat=2)
            binStatsExt.update(self.results[binId].geneCopyNumber(binIdToBinMarkerSets[binId]))
            fout.write(binId + '\t' + str(binStatsExt) + '\n')
        fout.close()

    def __writeMarkerGeneStats(self, directory, binIdToBinMarkerSets, bIndividualMarkers):
        markerGenesFile = os.path.join(directory, 'storage', DefaultValues.MARKER_GENE_STATS)
        fout = open(markerGenesFile, 'w')

        for binId in self.results:
            markerGenesHits = self.results[binId].getSummary(binIdToBinMarkerSets[binId], bIndividualMarkers, outputFormat=8)
            fout.write(binId + '\t' + str(markerGenesHits) + '\n')
        fout.close()

    def parseBinStats(self, resultsFolder, binStatsFile):
        """Read bin statistics from file."""
        binStatsFile = os.path.join(resultsFolder, 'storage', binStatsFile)

        checkFileExists(binStatsFile)

        binStats = {}
        with open(binStatsFile, 'r') as f:
            for line in f:
                lineSplit = line.split('\t')
                binStats[lineSplit[0]] = ast.literal_eval(lineSplit[1])

        return binStats

    def parseBinStatsExt(self, resultsFolder):
        """Read bin statistics from file."""
        binStatsExtFile = os.path.join(resultsFolder, 'storage', DefaultValues.BIN_STATS_EXT_OUT)

        checkFileExists(binStatsExtFile)

        binStatsExt = {}
        with open(binStatsExtFile, 'r') as f:
            for line in f:
                lineSplit = line.split('\t')
                binStatsExt[lineSplit[0]] = ast.literal_eval(lineSplit[1])

        return binStatsExt

    def parseMarkerGeneStats(self, resultsFolder):
        """Read bin statistics from file."""
        markerGeneStatsFile = os.path.join(resultsFolder, 'storage', DefaultValues.MARKER_GENE_STATS)

        checkFileExists(markerGeneStatsFile)

        markerGeneStats = {}
        with open(markerGeneStatsFile, 'r') as f:
            for line in f:
                lineSplit = line.split('\t')
                markerGeneStats[lineSplit[0]] = ast.literal_eval(lineSplit[1])

        return markerGeneStats

    def parseHmmerResults(self, fileName, resultsManager, bSkipAdjCorrection):
        """Parse HMMER results."""
        try:
            with open(fileName, 'r') as hmmerHandle:
                try:
                    HP = HMMERParser(hmmerHandle)
                except:
                    print("Error opening HMM file: ", fileName)
                    raise

                while True:
                    hit = HP.next()
                    if hit is None:
                        break
                    resultsManager.addHit(hit)

            # retain only best hit to PFAM clan
            pfam = PFAM(DefaultValues.PFAM_CLAN_FILE)
            resultsManager.markerHits = pfam.filterHitsFromSameClan(resultsManager.markerHits)

            # correct for errors in ORF calling
            if not bSkipAdjCorrection:
                resultsManager.identifyAdjacentMarkerGenes()

        except IOError as detail:
            sys.stderr.write(str(detail) + "\n")

    def __getHeader(self, outputFormat, binMarkerSets, coverageBinProfiles=None, table=None):
        """Get header for requested output table."""

        # keep count of single, double, triple genes etc...
        if outputFormat == 1:
            header = ['Bin Id', 'Marker lineage', '# genomes', '# markers', '# marker sets', '0', '1', '2', '3', '4', '5+', 'Completeness', 'Contamination', 'Strain heterogeneity']
        elif outputFormat == 2:
            header = ['Bin Id']
            header += ['Marker lineage', '# genomes', '# markers', '# marker sets']
            header += ['Completeness', 'Contamination', 'Strain heterogeneity']
            header += ['Genome size (bp)', '# ambiguous bases', '# scaffolds', '# contigs']
            header += ['N50 (scaffolds)', 'N50 (contigs)', 'Mean scaffold length (bp)', 'Mean contig length (bp)']
            header += ['Longest scaffold (bp)', 'Longest contig (bp)']
            header += ['GC', 'GC std (scaffolds > 1kbp)']
            header += ['Coding density', 'Translation table', '# predicted genes']
            header += ['0', '1', '2', '3', '4', '5+']

            if coverageBinProfiles != None:
                for bamId in coverageBinProfiles[list(coverageBinProfiles.keys())[0]]:
                    header += ['Coverage (' + bamId + ')', 'Coverage std (' + bamId + ')']

            if DefaultValues.MIN_SEQ_LEN_GC_STD != 1000:
                self.logger.error('Labeling error: GC std (scaffolds > 1kbp)')
                sys.exit(1)
        elif outputFormat == 3:
            header = ['Bin Id', 'Node Id', 'Marker lineage', '# genomes', '# markers', '# marker sets', '0', '1', '2', '3', '4', '5+', 'Completeness', 'Contamination', 'Strain heterogeneity']
        elif outputFormat == 4:
            header = None
        elif outputFormat == 5:
            header = ['Bin Id', 'Marker Id', 'Gene Id']
        elif outputFormat == 6:
            header = ['Bin Id', 'Marker Id', 'Gene Ids']
        elif outputFormat == 7:
            header = ['Bin Id', 'Marker Id', 'Gene Ids']
        elif outputFormat == 8:
            header = ['Bin Id', 'Gene Id', '{Marker Id, Start position, End position}']
        elif outputFormat == 9:
            if table is not None:
                header = ['Bin Id', 'Contig', 'Gene Number', 'Gene Start', 'Gene End', 'Gene Strand', 'Prot Length', 'Marker Id', 'Align Start', 'Align End', 'Sequence']
            else:
                header = " "
        elif outputFormat == 10:
            header = ['Scaffold Id', 'Bin Id', 'Length', '# contigs', 'GC', '# ORFs', 'Coding density', 'Marker Ids']

        return header

    def printSummary(self, outputFormat, aai, binIdToBinMarkerSets, bIndividualMarkers, coverageFile, bTabTable, outFile, anaFolder):
        # redirect output
        oldStdOut = reassignStdOut(outFile)

        coverageBinProfiles = None
        if coverageFile:
            coverage = Coverage(1)
            coverageBinProfiles = coverage.binProfiles(coverageFile)

        prettyTableFormats = [1, 2, 3, 9]

        header = self.__getHeader(outputFormat, binIdToBinMarkerSets[list(binIdToBinMarkerSets.keys())[0]], coverageBinProfiles, bTabTable)
        if bTabTable or outputFormat not in prettyTableFormats:
            bTabTable = True
            pTable = None

            if header != None:
                print('\t'.join(header))
        else:
            pTable = prettytable.PrettyTable(header)
            pTable.float_format = '.2'
            pTable.align = 'c'
            pTable.align[header[0]] = 'l'
            pTable.hrules = prettytable.FRAME
            pTable.vrules = prettytable.NONE

        seqsReported = 0
        for binId in sorted(self.results.keys()):
            seqsReported += self.results[binId].printSummary(outputFormat, aai, binIdToBinMarkerSets[binId], bIndividualMarkers, coverageBinProfiles, pTable, anaFolder)

        if outputFormat in [6, 7] and seqsReported == 0:
            print('[No marker genes satisfied the reporting criteria.]')

        if not bTabTable:
            if outputFormat in [1, 2]:
                print(pTable.get_string(sortby='Completeness', reversesort=True))
            else:
                # only print if there are rows
                if pTable.get_string(print_empty=False):
                    print(pTable.get_string(print_empty=False))

        # restore stdout
        restoreStdOut(outFile, oldStdOut) class HmmerAligner:
    def __init__(self, threads):
        self.logger = logging.getLogger('timestamp')
        self.totalThreads = threads

        self.outputFormat = 'Pfam'

    def makeAlignmentTopHit(self,
                               outDir,
                               hmmModelFile,
                               hmmTableFile,
                               binIdToModels,
                               bIgnoreThresholds,
                               evalueThreshold,
                               lengthThreshold,
                               bReportHitStats,
                               alignOutputDir,
                               bKeepUnmaskedAlign=False
                               ):
        """Align top hits in each bin. Assumes all bins are using the same marker genes."""

        self.logger.info("Extracting marker genes to align.")

        # parse HMM information
        resultsParser = ResultsParser(binIdToModels)

        # get HMM hits to each bin
        resultsParser.parseBinHits(outDir, hmmTableFile, False, bIgnoreThresholds, evalueThreshold, lengthThreshold)

        # extract the ORFs to align
        markerSeqs, markerStats = self.__extractMarkerSeqsTopHits(outDir, resultsParser)

        # generate individual HMMs required to create multiple sequence alignments
        binId = list(binIdToModels.keys())[0]
        hmmModelFiles = {}
        self.__makeAlignmentModels(hmmModelFile, binIdToModels[binId], hmmModelFiles)

        # align each of the marker genes
        makeSurePathExists(alignOutputDir)
        self.__alignMarkerGenes(markerSeqs, markerStats, bReportHitStats, hmmModelFiles, alignOutputDir, bKeepUnmaskedAlign)

        # remove the temporary HMM files
        for fileName in hmmModelFiles:
            os.remove(hmmModelFiles[fileName])

        return resultsParser

    def makeAlignmentToPhyloMarkers(self,
                                       outDir,
                                       hmmModelFile,
                                       hmmTableFile,
                                       binIdToModels,
                                       bIgnoreThresholds,
                                       evalueThreshold,
                                       lengthThreshold,
                                       bReportHitStats,
                                       alignOutputDir,
                                       bKeepUnmaskedAlign=False
                                       ):
        """Align hits to a set of common marker genes."""

        self.logger.info("Extracting marker genes to align.")

        # parse HMM information
        resultsParser = ResultsParser(binIdToModels)

        # get HMM hits to each bin
        resultsParser.parseBinHits(outDir, hmmTableFile, False, bIgnoreThresholds, evalueThreshold, lengthThreshold)

        # extract the ORFs to align
        markerSeqs, markerStats = self.__extractMarkerSeqsUnique(outDir, resultsParser)

        # generate individual HMMs required to create multiple sequence alignments
        binId = list(binIdToModels.keys())[0]
        hmmModelFiles = {}
        self.__makeAlignmentModels(hmmModelFile, binIdToModels[binId], hmmModelFiles)

        # align each of the marker genes
        makeSurePathExists(alignOutputDir)
        self.__alignMarkerGenes(markerSeqs, markerStats, bReportHitStats, hmmModelFiles, alignOutputDir, bKeepUnmaskedAlign)

        # remove the temporary HMM files
        for fileName in hmmModelFiles:
            os.remove(hmmModelFiles[fileName])

        return resultsParser

    def makeAlignmentsOfMultipleHits(self,
                                       outDir,
                                       markerFile,
                                       hmmTableFile,
                                       binIdToModels,
                                       binIdToBinMarkerSets,
                                       bIgnoreThresholds,
                                       evalueThreshold,
                                       lengthThreshold,
                                       alignOutputDir,
                                       ):
        """Align markers with multiple hits within a bin."""

        makeSurePathExists(alignOutputDir)

        # parse HMM information
        resultsParser = ResultsParser(binIdToModels)

        # get HMM hits to each bin
        resultsParser.parseBinHits(outDir, hmmTableFile, False, bIgnoreThresholds, evalueThreshold, lengthThreshold)

        # align any markers with multiple hits in a bin
        self.logger.info('Aligning marker genes with multiple hits in a single bin:')

        # process each bin in parallel
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        for binId in binIdToModels:
            workerQueue.put(binId)

        for _ in range(self.totalThreads):
            workerQueue.put(None)

        try:
            calcProc = [mp.Process(target=self.__createMSA, args=(resultsParser, binIdToBinMarkerSets, markerFile, outDir, alignOutputDir, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__reportBinProgress, args=(len(binIdToModels), writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put(None)
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

    def __createMSA(self, resultsParser, binIdToBinMarkerSets, hmmModelFile, outDir, alignOutputDir, queueIn, queueOut):
        """Create multiple sequence alignment for markers with multiple hits in a bin."""

        HF = HMMERRunner(mode='fetch')

        while True:
            binId = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            markersWithMultipleHits = self.__extractMarkersWithMultipleHits(outDir, binId, resultsParser, binIdToBinMarkerSets[binId])

            if len(markersWithMultipleHits) != 0:
                # create multiple sequence alignments for markers with multiple hits
                binAlignOutputDir = os.path.join(alignOutputDir, binId)
                makeSurePathExists(binAlignOutputDir)
                for markerId in markersWithMultipleHits:
                    tempModelFile = os.path.join(tempfile.gettempdir(), str(uuid.uuid4()))
                    HF.fetch(hmmModelFile, markerId, tempModelFile)

                    self.__alignMarker(markerId, markersWithMultipleHits[markerId], None, False, binAlignOutputDir, tempModelFile, bKeepUnmaskedAlign=False)

                    os.remove(tempModelFile)

            queueOut.put(binId)

    def __reportBinProgress(self, numBins, queueIn):
        """Report number of processed bins."""

        numProcessedBins = 0
        if self.logger.getEffectiveLevel() <= logging.INFO:
            statusStr = '    Finished processing %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
            sys.stderr.write('%s\r' % statusStr)
            sys.stderr.flush()

        while True:
            binId = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            if self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedBins += 1
                statusStr = '    Finished processing %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

    def __alignMarkerGenes(self, markerSeqs, markerStats, bReportHitStats, hmmModelFiles, alignOutputDir, bKeepUnmaskedAlign=False, bReportProgress=True):
        """Align marker genes with HMMs in parallel."""

        if bReportProgress:
            self.logger.info("Aligning %d marker genes with %d threads:" % (len(hmmModelFiles), self.totalThreads))

        # process each bin in parallel
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        for markerId in hmmModelFiles:
            workerQueue.put(markerId)

        for _ in range(self.totalThreads):
            workerQueue.put(None)

        try:
            calcProc = [mp.Process(target=self.__alignMarkerParallel, args=(markerSeqs, markerStats, bReportHitStats, alignOutputDir, hmmModelFiles, bKeepUnmaskedAlign, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__reportAlignmentProgress, args=(len(hmmModelFiles), bReportProgress, writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put(None)
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

    def __alignMarkerParallel(self, markerSeqs, markerStats, bReportHitStats, alignOutputDir, hmmModelFiles, bKeepUnmaskedAlign, queueIn, queueOut):
        while True:
            markerId = queueIn.get(block=True, timeout=None)
            if markerId == None:
                break

            self.__alignMarker(markerId, markerSeqs[markerId], markerStats[markerId], bReportHitStats, alignOutputDir, hmmModelFiles[markerId], bKeepUnmaskedAlign)

            queueOut.put(markerId)

    def __alignMarker(self, markerId, binSeqs, binStats, bReportHitStats, alignOutputDir, hmmModelFile, bKeepUnmaskedAlign):
        unalignSeqFile = os.path.join(alignOutputDir, markerId + '.unaligned.faa')
        fout = open(unalignSeqFile, 'w')
        numSeqs = 0
        for binId, seqs in binSeqs.items():
            for seqId, seq in seqs.items():
                header = '>' + binId + DefaultValues.SEQ_CONCAT_CHAR + seqId
                if bReportHitStats:
                    header += ' [e-value=%.4g,score=%.1f]' % (binStats[binId][seqId][0], binStats[binId][seqId][1])

                fout.write(header + '\n')
                fout.write(seq + '\n')
                numSeqs += 1
        fout.close()

        if numSeqs > 0:
            alignSeqFile = os.path.join(alignOutputDir, markerId + '.aligned.faa')
            HA = HMMERRunner(mode='align')
            HA.align(hmmModelFile, unalignSeqFile, alignSeqFile, writeMode='>', outputFormat=self.outputFormat, trim=False)

            makedSeqFile = os.path.join(alignOutputDir, markerId + '.masked.faa')
            self.__maskAlignment(alignSeqFile, makedSeqFile)

            if not bKeepUnmaskedAlign:
                os.remove(alignSeqFile)

        os.remove(unalignSeqFile)

    def __reportAlignmentProgress(self, numMarkers, bReportProgress, queueIn):
        """Report number of processed markers."""

        numProcessedGenes = 0
        if bReportProgress and self.logger.getEffectiveLevel() <= logging.INFO:
            statusStr = '    Finished aligning %d of %d (%.2f%%) marker genes.' % (numProcessedGenes, numMarkers, float(numProcessedGenes) * 100 / numMarkers)
            sys.stderr.write('%s\r' % statusStr)
            sys.stderr.flush()

        while True:
            binId = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            if bReportProgress and self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedGenes += 1
                statusStr = '    Finished aligning %d of %d (%.2f%%) marker genes.' % (numProcessedGenes, numMarkers, float(numProcessedGenes) * 100 / numMarkers)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        if bReportProgress and self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

    def __maskAlignment(self, inputFile, outputFile):
        """Read HMMER alignment in STOCKHOLM format and output masked alignment in FASTA format."""
        # read STOCKHOLM alignment
        seqs = {}
        seqStats = {}
        for line in open(inputFile):
            line = line.rstrip()
            if line == '' or line[0] == '#' or line == '//':
                if 'GC RF' in line:
                    mask = line.split('GC RF')[1].strip()
                elif '=GS' in line:
                    # read additional sequence informations
                    lineSplit = line.split()
                    seqId = lineSplit[1]
                    stats = lineSplit[3].strip()
                    seqStats[seqId] = stats
                continue
            else:
                lineSplit = line.split()
                seqs[lineSplit[0]] = lineSplit[1].upper().replace('.', '-').strip()

        # output masked sequences in FASTA format
        fout = open(outputFile, 'w')
        for seqId, seq in seqs.items():
            if seqStats:
                fout.write('>%s %s\n' % (seqId, seqStats[seqId]))
            else:
                fout.write('>' + seqId + '\n')

            maskedSeq = ''.join([seq[i] for i in range(0, len(seq)) if mask[i] == 'x'])
            fout.write(maskedSeq + '\n')
        fout.close()

    def __extractMarkerSeqsTopHits(self, outDir, resultsParser):
        """Extract marker sequences from top hits (assume all bins use the same HMM file)."""

        markerSeqs = defaultdict(dict)
        markerStats = defaultdict(dict)
        for binId in resultsParser.results:
            # read ORFs for bin
            aaGeneFile = os.path.join(outDir, 'bins', binId, DefaultValues.PRODIGAL_AA)
            binORFs = readFasta(aaGeneFile)

            # extract ORFs hitting a marker
            for markerId, hits in resultsParser.results[binId].markerHits.items():
                markerSeqs[markerId][binId] = {}
                markerStats[markerId][binId] = {}

                # sort hits from highest to lowest e-value in order to ensure only the best hit
                # to a given target is retained
                hits.sort(key=lambda x: x.full_e_value, reverse=True)
                topHit = hits[0]
                markerSeqs[markerId][binId][topHit.target_name] = self.__extractSeq(topHit.target_name, binORFs)
                markerStats[markerId][binId][topHit.target_name] = [topHit.full_e_value, topHit.full_score]

        return markerSeqs, markerStats

    def __extractMarkerSeqsUnique(self, outDir, resultsParser):
        """Extract marker sequences with a single unique hit."""

        markerSeqs = defaultdict(dict)
        markerStats = defaultdict(dict)
        for binId in resultsParser.results:
            # read ORFs for bin
            aaGeneFile = os.path.join(outDir, 'bins', binId, DefaultValues.PRODIGAL_AA)
            binORFs = readFasta(aaGeneFile)

            # extract ORFs hitting a marker
            for markerId, hits in resultsParser.results[binId].markerHits.items():
                markerSeqs[markerId][binId] = {}
                markerStats[markerId][binId] = {}

                # only record hits which are unique
                if len(hits) == 1:
                    hit = hits[0]
                    markerSeqs[markerId][binId][hit.target_name] = self.__extractSeq(hit.target_name, binORFs)
                    markerStats[markerId][binId][hit.target_name] = [hit.full_e_value, hit.full_score]

        return markerSeqs, markerStats

    def __extractSeq(self, seqId, seqs):
        """Extract sequence data."""

        if DefaultValues.SEQ_CONCAT_CHAR in seqId:
            seqIds = seqId.split(DefaultValues.SEQ_CONCAT_CHAR)

            seq = ''
            for seqId in seqIds:
                tempSeq = seqs[seqId]
                if tempSeq[-1] == '*':
                    tempSeq = tempSeq[0:-1]  # remove final '*' inserted by prodigal

                seq += tempSeq

            rtnSeq = seq
        else:
            rtnSeq = seqs[seqId]

            if rtnSeq[-1] == '*':
                rtnSeq = rtnSeq[0:-1]  # remove final '*' inserted by prodigal

        return rtnSeq

    def __extractMarkersWithMultipleHits(self, outDir, binId, resultsParser, binMarkerSet):
        """Extract markers with multiple hits within a single bin."""

        markersWithMultipleHits = defaultdict(dict)

        aaGeneFile = os.path.join(outDir, 'bins', binId, DefaultValues.PRODIGAL_AA)
        binORFs = readFasta(aaGeneFile)

        markerGenes = binMarkerSet.selectedMarkerSet().getMarkerGenes()
        for markerId, hits in resultsParser.results[binId].markerHits.items():
            if markerId not in markerGenes or len(hits) < 2:
                continue

            # sort hits from highest to lowest e-value in order to ensure only the best hit
            # to a given target is retained
            hits.sort(key=lambda x: x.full_e_value, reverse=True)

            # Note: this data structure is used to mimic that used by __extractMarkerSeqsTopHits()
            markersWithMultipleHits[markerId][binId] = {}
            for hit in hits:
                markersWithMultipleHits[markerId][binId][hit.target_name] = self.__extractSeq(hit.target_name, binORFs)

        return markersWithMultipleHits

    def __makeAlignmentModels(self, hmmModelFile, modelKeys, hmmModelFiles, bReportProgress=True):
        """Make temporary HMM files used to create HMM alignments."""

        if bReportProgress:
            self.logger.info("Extracting %d HMMs with %d threads:" % (len(modelKeys), self.totalThreads))

        # process each marker in parallel
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        for modelId in modelKeys:
            fetchFilename = os.path.join(tempfile.gettempdir(), str(uuid.uuid4()))
            hmmModelFiles[modelId] = fetchFilename
            workerQueue.put((modelId, fetchFilename))

        for _ in range(self.totalThreads):
            workerQueue.put((None, None))

        try:
            calcProc = [mp.Process(target=self.__extractModel, args=(hmmModelFile, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__reportModelExtractionProgress, args=(len(modelKeys), bReportProgress, writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put(None)
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

    def __extractModel(self, hmmModelFile, queueIn, queueOut):
        """Extract HMM."""
        HF = HMMERRunner(mode='fetch')

        while True:
            modelId, fetchFilename = queueIn.get(block=True, timeout=None)
            if modelId == None:
                break

            HF.fetch(hmmModelFile, modelId, fetchFilename)

            queueOut.put(modelId)

    def __reportModelExtractionProgress(self, numMarkers, bReportProgress, queueIn):
        """Report number of extracted HMMs."""

        numModelsExtracted = 0
        if bReportProgress and self.logger.getEffectiveLevel() <= logging.INFO:
            statusStr = '    Finished extracting %d of %d (%.2f%%) HMMs.' % (numModelsExtracted, numMarkers, float(numModelsExtracted) * 100 / numMarkers)
            sys.stderr.write('%s\r' % statusStr)
            sys.stderr.flush()

        while True:
            modelId = queueIn.get(block=True, timeout=None)
            if modelId == None:
                break

            if bReportProgress and self.logger.getEffectiveLevel() <= logging.INFO:
                numModelsExtracted += 1
                statusStr = '    Finished extracting %d of %d (%.2f%%) HMMs.' % (numModelsExtracted, numMarkers, float(numModelsExtracted) * 100 / numMarkers)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        if bReportProgress and self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')
 class MarkerGeneFinder():
    """Identify marker genes within binned sequences using Prodigal and HMMER."""
    def __init__(self, threads):
        self.logger = logging.getLogger('timestamp')
        self.totalThreads = threads

    def find(self, binFiles, outDir, tableOut, hmmerOut, markerFile, bKeepAlignment, bNucORFs, bCalledGenes):
        """Identify marker genes in each bin using prodigal and HMMER."""

        # make sure HMMER and prodigal are on system path
        HMMERRunner()

        if not bCalledGenes:
            ProdigalRunner('')

        # process each fasta file
        self.threadsPerSearch = max(1, int(self.totalThreads / len(binFiles)))
        self.logger.info("Identifying marker genes in %d bins with %d threads:" % (len(binFiles), self.totalThreads))

        # process each bin in parallel
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        for binFile in binFiles:
            workerQueue.put(binFile)

        for _ in range(self.totalThreads):
            workerQueue.put(None)

        binIdToModels = mp.Manager().dict()

        try:
            calcProc = [mp.Process(target=self.__processBin, args=(outDir, tableOut, hmmerOut, markerFile, bKeepAlignment, bNucORFs, bCalledGenes, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__reportProgress, args=(len(binFiles), binIdToModels, writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put((None, None))
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

        # create a standard dictionary from the managed dictionary
        d = {}
        for binId in binIdToModels.keys():
            d[binId] = binIdToModels[binId]

        return d

    def __processBin(self, outDir, tableOut, hmmerOut, markerFile, bKeepAlignment, bNucORFs, bCalledGenes, queueIn, queueOut):
        """Thread safe bin processing."""

        markerSetParser = MarkerSetParser(self.threadsPerSearch)

        while True:
            binFile = queueIn.get(block=True, timeout=None)
            if binFile == None:
                break

            binId = binIdFromFilename(binFile)
            binDir = os.path.join(outDir, 'bins', binId)
            makeSurePathExists(binDir)

            # run Prodigal
            if not bCalledGenes:
                prodigal = ProdigalRunner(binDir)
                if not prodigal.areORFsCalled(bNucORFs):
                    prodigal.run(binFile, bNucORFs)
                aaGeneFile = prodigal.aaGeneFile
            else:
                aaGeneFile = binFile
                shutil.copyfile(aaGeneFile, os.path.join(binDir, DefaultValues.PRODIGAL_AA))

            # extract HMMs into temporary file
            hmmModelFile = markerSetParser.createHmmModelFile(binId, markerFile)

            # run HMMER
            hmmer = HMMERRunner()
            tableOutPath = os.path.join(binDir, tableOut)
            hmmerOutPath = os.path.join(binDir, hmmerOut)

            keepAlignStr = ''
            if not bKeepAlignment:
                keepAlignStr = '--noali'
            hmmer.search(hmmModelFile, aaGeneFile, tableOutPath, hmmerOutPath,
                         '--cpu ' + str(self.threadsPerSearch) + ' --notextw -E 0.1 --domE 0.1 ' + keepAlignStr,
                         bKeepAlignment)

            queueOut.put((binId, hmmModelFile))

    def __reportProgress(self, numBins, binIdToModels, queueIn):
        """Report number of processed bins."""

        numProcessedBins = 0
        if self.logger.getEffectiveLevel() <= logging.INFO:
            statusStr = '    Finished processing %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
            sys.stderr.write('%s\r' % statusStr)
            sys.stderr.flush()

        while True:
            binId, hmmModelFile = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            # parse HMM file
            # (This is done here as pushing the models onto the shared queue is too memory intensive)
            modelParser = HmmModelParser(hmmModelFile)
            models = modelParser.models()

            binIdToModels[binId] = models

            if os.path.exists(hmmModelFile):
                os.remove(hmmModelFile)

            indexFile = hmmModelFile + '.ssi'
            if os.path.exists(indexFile):
                os.remove(indexFile)

            if self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedBins += 1
                statusStr = '    Finished processing %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n') class PplacerRunner():
    """Wrapper for running pplacer."""
    def __init__(self, threads):
        self.logger = logging.getLogger('timestamp')
        self.numThreads = threads

        # make sure pplace and guppy are on the system path
        self.__checkForPplacer()
        self.__checkForGuppy()

    def run(self, binFiles, resultsParser, outDir, bReducedTree):
        # make sure output and tree directories exist
        checkDirExists(outDir)
        alignOutputDir = os.path.join(outDir, 'storage', 'tree')
        checkDirExists(alignOutputDir)

        treeFile = os.path.join(alignOutputDir, DefaultValues.PPLACER_TREE_OUT)
        pplacerJsonOut = os.path.join(alignOutputDir, DefaultValues.PPLACER_JSON_OUT)
        pplacerOut = os.path.join(alignOutputDir, DefaultValues.PPLACER_OUT)

        # create concatenated alignment file for each bin
        concatenatedAlignFile = self.__createConcatenatedAlignment(binFiles, resultsParser, alignOutputDir)

        pplacerRefPkg = DefaultValues.PPLACER_REF_PACKAGE_FULL
        if bReducedTree:
            pplacerRefPkg = DefaultValues.PPLACER_REF_PACKAGE_REDUCED

        # check if concatenated alignment file is empty
        # (this can occur when all genomes have no phylogenetically informative marker genes)
        if os.stat(concatenatedAlignFile)[stat.ST_SIZE] == 0:
            self.logger.info('No genomes were identified that could be placed in the reference genome tree.')
            shutil.copyfile(os.path.join(pplacerRefPkg, DefaultValues.GENOME_TREE), treeFile)
            return

        # run pplacer to place bins in reference genome tree
        self.logger.info('Placing %d bins into the genome tree with pplacer (be patient).' % len(binFiles))
        cmd = 'pplacer -j %d -c %s -o %s %s > %s' % (self.numThreads,
                                                     pplacerRefPkg,
                                                     pplacerJsonOut,
                                                     concatenatedAlignFile,
                                                     pplacerOut)
        os.system(cmd)

        # extract tree
        cmd = 'guppy tog -o %s %s' % (treeFile, pplacerJsonOut)
        os.system(cmd)

    def __createConcatenatedAlignment(self, binFiles, resultsParser, alignOutputDir):
        """Create a concatenated alignment of marker genes for each bin."""

        # read alignment files
        self.logger.info('Reading marker alignment files.')
        alignments = defaultdict(dict)
        files = os.listdir(alignOutputDir)
        binIds = set()
        for f in files:
            if f.endswith('.masked.faa'):
                markerId = f[0:f.find('.masked.faa')]
                seqs = readFasta(os.path.join(alignOutputDir, f))

                for seqId, seq in seqs.items():
                    binId = seqId[0:seqId.find(DefaultValues.SEQ_CONCAT_CHAR)]

                    alignments[markerId][binId] = seq
                    binIds.add(binId)

        # get all markers and their lengths
        markerIds = resultsParser.models[list(resultsParser.models.keys())[0]].keys()
        markerIdLens = {}
        for markerId in markerIds:
            markerIdLens[markerId] = resultsParser.models[list(resultsParser.models.keys())[0]][markerId].leng

        # create concatenated alignment
        self.logger.info('Concatenating alignments.')
        concatenatedSeqs = {}
        for markerId in sorted(markerIds):
            seqs = alignments[markerId]

            for binId in binIds:
                if binId in seqs:
                    # append alignment
                    concatenatedSeqs[binId] = concatenatedSeqs.get(binId, '') + seqs[binId]
                else:
                    # missing gene
                    concatenatedSeqs[binId] = concatenatedSeqs.get(binId, '') + '-' * markerIdLens[markerId]

        # save concatenated alignment
        concatenatedAlignFile = os.path.join(alignOutputDir, DefaultValues.PPLACER_CONCAT_SEQ_OUT)
        writeFasta(concatenatedSeqs, concatenatedAlignFile)

        return concatenatedAlignFile

    def __checkForPplacer(self):
        """Check to see if pplacer is on the system before we try to run it."""

        # Assume that a successful pplacer -h returns 0 and anything
        # else returns something non-zero
        try:
            subprocess.call(['pplacer', '-h'], stdout=open(os.devnull, 'w'), stderr=subprocess.STDOUT)
        except:
            self.logger.error("Make sure pplacer is on your system path.")
            sys.exit(1)

    def __checkForGuppy(self):
        """Check to see if guppy is on the system before we try to run it."""

        # Assume that a successful pplacer -h returns 0 and anything
        # else returns something non-zero
        try:
            subprocess.call(['guppy', '-h'], stdout=open(os.devnull, 'w'), stderr=subprocess.STDOUT)
        except:
            self.logger.error("Make sure guppy, which is part of the pplacer package, is on your system path.")
            sys.exit(1)
 class TreeParser():
    """Parse genome tree and associated tree metadata."""
    def __init__(self):
        self.logger = logging.getLogger('timestamp')

    def printSummary(self, outputFormat, outDir, resultsParser, bTabTable, outFile, binStats):
        if outputFormat == 1:
            self.reportBinTaxonomy(outDir, resultsParser, bTabTable, outFile, binStats, bLineageStatistics=False)
        elif outputFormat == 2:
            self.reportBinTaxonomy(outDir, resultsParser, bTabTable, outFile, binStats, bLineageStatistics=True)
        elif outputFormat == 3:
            self.reportNewickTree(outDir, outFile, None)
        elif outputFormat == 4:
            self.reportNewickTree(outDir, outFile, 'taxonomy')
        elif outputFormat == 5:
            self.reportFullMSA(outDir, outFile)
        else:
            self.logger.error("Unknown output format: %d", outputFormat)

    def reportFullMSA(self, outDir, outFile):
        """Create MSA with all reference and bin alignments."""

        # write bin alignments to file
        oldStdOut = reassignStdOut(outFile)
        for line in open(os.path.join(outDir, 'storage', 'tree', DefaultValues.PPLACER_CONCAT_SEQ_OUT)):
            print(line.rstrip())

        # read duplicate seqs
        duplicateNodes = self.__readDuplicateSeqs()

        # write reference alignments to file
        seqs = readFasta(os.path.join(DefaultValues.PPLACER_REF_PACKAGE_FULL, DefaultValues.GENOME_TREE_FASTA))
        for seqId, seq in seqs.items():
            print('>' + seqId)
            print(seq)

            if seqId in duplicateNodes:
                for dupSeqId in duplicateNodes[seqId]:
                    print('>' + dupSeqId)
                    print(seq)

        restoreStdOut(outFile, oldStdOut)

    def readPlacementFile(self, placementFile):
        '''Read pplacer JSON placement file.'''
        jsonData = open(placementFile)

        data = json.load(jsonData)
        binIdToPP = {}
        for placementData in data['placements']:
            binId = placementData['nm'][0][0]

            topPP = 0
            for pp in placementData['p']:
                if pp[2] > topPP:
                    topPP = pp[2]

            binIdToPP[binId] = topPP

        jsonData.close()

        return binIdToPP

    def reportNewickTree(self, outDir, outFile, leafLabels=None):
        # read duplicate nodes
        duplicateSeqs = self.__readDuplicateSeqs()

        # read tree
        treeFile = os.path.join(outDir, 'storage', 'tree', DefaultValues.PPLACER_TREE_OUT)
        tree = dendropy.Tree.get_from_path(treeFile, schema='newick', rooting="force-rooted", preserve_underscores=True)

        # clean up internal node labels
        for node in tree.internal_nodes():
            if node.label:
                labelSplit = node.label.split('|')

                label = labelSplit[0]
                if labelSplit[1] != '':
                    label += '|' + labelSplit[1]
                if labelSplit[2] != '':
                    label += '|' + labelSplit[2]

                node.label = label

        # insert duplicate nodes into tree
        for leaf in tree.leaf_nodes():
            duplicates = duplicateSeqs.get(leaf.taxon.label, None)
            if duplicates != None:
                newParent = leaf.parent_node.new_child(edge_length=leaf.edge_length)
                curLeaf = leaf.parent_node.remove_child(leaf)
                newParent.new_child(taxon=curLeaf.taxon, edge_length=0)
                for d in duplicates:
                    newParent.new_child(taxon=dendropy.Taxon(label=d), edge_length=0)

        # append taxonomy to leaf nodes
        if leafLabels == 'taxonomy':
            # read taxonomy string for each IMG genome
            taxonomy = {}
            for line in open(os.path.join(DefaultValues.GENOME_TREE_DIR, DefaultValues.GENOME_TREE_TAXONOMY)):
                lineSplit = line.split('\t')
                taxonomy[lineSplit[0]] = lineSplit[1].rstrip()

            # append taxonomy to leaf labels
            for leaf in tree.leaf_nodes():
                taxaStr = taxonomy.get(leaf.taxon.label, None)
                if taxaStr:
                    leaf.taxon.label += '|' + taxaStr

        # write out tree
        oldStdOut = reassignStdOut(outFile)
        print(tree.as_string(schema='newick', suppress_rooting=True))
        restoreStdOut(outFile, oldStdOut)

    def getInsertionBranchId(self, outDir, binIds):
        # make sure output and tree directories exist
        checkDirExists(outDir)
        alignOutputDir = os.path.join(outDir, 'storage', 'tree')
        checkDirExists(alignOutputDir)

        # read genome tree (if it exists)
        binIdToUID = {}
        treeFile = os.path.join(alignOutputDir, DefaultValues.PPLACER_TREE_OUT)
        tree = dendropy.Tree.get_from_path(treeFile, schema='newick', rooting="force-rooted", preserve_underscores=True)

        # find first parent of each bin with a taxonomic label
        for binId in binIds:
            node = tree.find_node_with_taxon_label(binId)
            if node == None:
                binIdToUID[binId] = 'NA'
                continue

            # find first node decorated with a UID string between leaf and root
            parentNode = node.parent_node
            while parentNode != None:
                if parentNode.label:
                    uid = parentNode.label.split('|')[0]
                    break

                parentNode = parentNode.parent_node

            binIdToUID[binId] = uid

        return binIdToUID

    def getBinTaxonomy(self, outDir, binIds):
        # make sure output and tree directories exist
        checkDirExists(outDir)
        alignOutputDir = os.path.join(outDir, 'storage', 'tree')
        checkDirExists(alignOutputDir)

        # read genome tree (if it exists)
        binIdToTaxonomy = {}
        treeFile = os.path.join(alignOutputDir, DefaultValues.PPLACER_TREE_OUT)
        tree = dendropy.Tree.get_from_path(treeFile, schema='newick', rooting="force-rooted", preserve_underscores=True)

        # find first parent of each bin with a taxonomic label
        for binId in binIds:
            node = tree.find_node_with_taxon_label(binId)
            if node == None:
                binIdToTaxonomy[binId] = 'NA'
                continue

            # find first node decorated with a taxon string between leaf and root
            taxaStr = None
            parentNode = node.parent_node
            while parentNode != None:
                if parentNode.label:
                    tokens = parentNode.label.split('|')

                    if tokens[1] != '':
                        if taxaStr:
                            taxaStr = tokens[1] + ';' + taxaStr
                        else:
                            taxaStr = tokens[1]

                parentNode = parentNode.parent_node

            if not taxaStr:
                domainNode = self.__findDomainNode(node)
                taxaStr = domainNode.label.split('|')[1] + ' (root)'

            binIdToTaxonomy[node.taxon.label] = taxaStr

        return binIdToTaxonomy

    def __findDomainNode(self, binNode):
        """Find node defining the domain. Assumes 'binNode' is the leaf node of a bin on either the archaeal or bacterial branch."""

        # bin is either on the bacterial or archaeal branch descendant from the root,
        # so descend tree to first internal node with a label. Note that there may
        # be multiple bins inserted in the tree creating unlabeled internal nodes.
        # Identical population bins can also cause internal bifurcations.

        # find first parent node which contains at least one IMG genome
        curNode = binNode.parent_node
        while True:
            found_ref_genome = False
            for leaf in curNode.leaf_nodes():
                if leaf.taxon.label.startswith('IMG_'):
                    found_ref_genome = True
                    break
                    
            if found_ref_genome:
                break

            curNode = curNode.parent_node

        # perform depth first search to find a labeled internal node
        queue = [curNode]
        while queue:
            curNode = queue.pop(0)

            if curNode.label:
                return curNode

            for child in curNode.child_nodes():
                if child.is_internal():
                    queue.append(child)

        self.logger.error('Failed to associate bin with a domain. Please report this bug.')
        sys.exit(1)

    def getBinSisterTaxonomy(self, outDir, binIds):
        # make sure output and tree directories exist
        checkDirExists(outDir)
        alignOutputDir = os.path.join(outDir, 'storage', 'tree')
        checkDirExists(alignOutputDir)

        # read genome tree
        treeFile = os.path.join(alignOutputDir, DefaultValues.PPLACER_TREE_OUT)
        tree = dendropy.Tree.get_from_path(treeFile, schema='newick', rooting="force-rooted", preserve_underscores=True)

        # read taxonomy string for each IMG genome
        leafIdToTaxonomy = {}
        for line in open(os.path.join(DefaultValues.GENOME_TREE_DIR, DefaultValues.GENOME_TREE_TAXONOMY)):
            lineSplit = line.split('\t')
            leafIdToTaxonomy[lineSplit[0]] = lineSplit[1].rstrip()

        # find LCA of all labeled node in sister lineage
        binIdToSisterTaxonomy = {}
        for binId in binIds:
            node = tree.find_node_with_taxon_label(binId)

            taxaStr = ''
            if node != None:
                # get taxonomic labels of all internal nodes in sister lineages
                sisterNodes = node.sister_nodes()
                internalTaxonomyLabels = set()
                leafTaxonomyLabels = set()
                for sn in sisterNodes:
                    for curNode in sn.postorder_iter():
                        if curNode.is_leaf():
                            if curNode.taxon.label:
                                taxonomy = leafIdToTaxonomy.get(curNode.taxon.label, None)
                                if taxonomy != None:  # inserted bins will not have an assigned taxonomy
                                    for taxa in taxonomy.split(';'):
                                        leafTaxonomyLabels.add(taxa.strip())
                        else:
                            if curNode.label:
                                tokens = curNode.label.split('|')
                                if tokens[1] != '':
                                    for taxa in tokens[1].split(';'):
                                        internalTaxonomyLabels.add(taxa)

                # find LCA of taxonomic labels in rank order;
                # only consider leaf node labels if there were no internal labels
                labels = internalTaxonomyLabels
                if len(labels) == 0:
                    labels = leafTaxonomyLabels

                for prefix in taxonomicPrefixes:
                    taxa = []
                    for taxon in labels:
                        if prefix in taxon:
                            taxa.append(taxon)

                    if len(taxa) == 1:
                        # unambiguous label at this rank
                        taxaStr += taxa[0] + ';'
                    elif len(taxa) > 1:
                        # unable to resolve taxonomy at this rank
                        break

            if not taxaStr:
                taxaStr = 'unresolved'
            binIdToSisterTaxonomy[binId] = taxaStr

        return binIdToSisterTaxonomy

    def __getNextNamedNode(self, node, uniqueIdToLineageStatistics):
        """Get first parent node with taxonomy information."""
        parentNode = node.parent_node
        while True:
            if parentNode == None:
                break  # reached the root node so terminate

            if parentNode.label:
                trustedUniqueId = parentNode.label.split('|')[0]
                trustedStats = uniqueIdToLineageStatistics[trustedUniqueId]
                if trustedStats['taxonomy'] != '':
                    return trustedStats['taxonomy']

            parentNode = parentNode.parent_node

        return 'root'

    def __getMarkerSet(self, parentNode, tree, uniqueIdToLineageStatistics,
                                    numGenomesMarkers, bootstrap,
                                    bForceDomain, bRequireTaxonomy):
        """Get marker set for next parent node meeting selection criteria."""

        # ascend tree to root finding first node meeting all selection criteria
        selectedParentNode = parentNode
        taxonomyStr = 'root'
        while True:
            if selectedParentNode.label:  # nodes inserted by PPLACER will not have a label
                trustedUniqueId = selectedParentNode.label.split('|')[0]
                nodeTaxonomy = selectedParentNode.label.split('|')[1]

                stats = uniqueIdToLineageStatistics[trustedUniqueId]
                if ((stats['# genomes'] == 'NA' 
                    or int(stats['# genomes']) >= numGenomesMarkers) 
                    and (stats['bootstrap'] == 'NA' or int(stats['bootstrap']) >= bootstrap)):
                    if not bForceDomain or nodeTaxonomy in ['k__Bacteria', 'k__Archaea']:
                        if not bRequireTaxonomy or stats['taxonomy'] != '':
                            # get closest taxonomic label
                            taxonomyStr = stats['taxonomy']
                            if not bRequireTaxonomy and stats['taxonomy'] == '':
                                taxonomyStr = self.__getNextNamedNode(selectedParentNode, uniqueIdToLineageStatistics)

                            # all criteria meet, so use marker set from this node
                            break

            if selectedParentNode.parent_node == None:
                break  # reached the root node so terminate

            selectedParentNode = selectedParentNode.parent_node

        # get marker set meeting all criteria required for a trusted marker set
        taxonomyStr = taxonomyStr.split(';')[-1]  # extract most specific taxonomy identifier
        markerSet = MarkerSet(trustedUniqueId, taxonomyStr, int(stats['# genomes']), eval(stats['marker set']))

        return selectedParentNode, markerSet

    def __refineMarkerSet(self, markerSet, binNode, tree, uniqueIdToLineageStatistics, numGenomesRefine):
        """Refine marker set to account for lineage-specific gene loss and duplication."""

        # lineage-specific refine is done with the sister lineage to where the bin is inserted

        # get lineage-specific marker set which will be used to refine the above marker set
        curNode = binNode.sister_nodes()[0]
        while True:
            if curNode.label:  # nodes inserted by PPLACER will not have a label
                uniqueId = curNode.label.split('|')[0]
                stats = uniqueIdToLineageStatistics[uniqueId]

                if stats['# genomes'] >= numGenomesRefine:
                    break

            curNode = curNode.parent_node
            if curNode == None:
                break  # reached the root node so terminate

        # get lineage-specific marker set
        lineageMarkerSet = eval(stats['marker set'])

        # refine marker set by finding the intersection between these two sets,
        # this removes markers that are not single-copy or ubiquitous in the
        # specific lineage of a bin
        # Note: co-localization information is taken from the trusted set

        # get all lineage-specific marker genes
        allLineageSpecificGenes = set()
        for m in lineageMarkerSet:
            for gene in m:
                allLineageSpecificGenes.add(gene)

        # remove genes not present in the lineage-specific gene set
        finalMarkerSet = []
        for ms in markerSet.markerSet:
            s = set()
            for gene in ms:
                if gene in allLineageSpecificGenes:
                    s.add(gene)

            if s:
                finalMarkerSet.append(s)

        refinedMarkerSet = MarkerSet(markerSet.UID, markerSet.lineageStr, markerSet.numGenomes, finalMarkerSet)

        return refinedMarkerSet

    def __readLineageSpecificGenesToRemove(self):
        """Get set of genes subject to lineage-specific gene loss and duplication."""

        self.lineageSpecificGenesToRemove = {}
        for line in open(os.path.join(DefaultValues.GENOME_TREE_DIR, DefaultValues.GENOME_TREE_MISSING_DUPLICATE)):
            lineSplit = line.split('\t')
            uid = lineSplit[0]
            missingGenes = eval(lineSplit[1])
            duplicateGenes = eval(lineSplit[2])
            self.lineageSpecificGenesToRemove[uid] = missingGenes.union(duplicateGenes)

    def __removeInvalidLineageMarkerGenes(self, markerSet, lineageSpecificMarkersToRemove):
        """Refine marker set to account for lineage-specific gene loss and duplication."""

        # refine marker set by removing marker genes subject to lineage-specific
        # gene loss and duplication
        #
        # Note: co-localization information is taken from the trusted set

        finalMarkerSet = []
        for ms in markerSet.markerSet:
            s = set()
            for gene in ms:
                geneIdToTest = gene
                if geneIdToTest.startswith('PF'):
                    geneIdToTest = gene.replace('PF', 'pfam')
                    geneIdToTest = geneIdToTest[0:geneIdToTest.rfind('.')]

                if geneIdToTest not in lineageSpecificMarkersToRemove:
                    s.add(gene)

            if s:
                finalMarkerSet.append(s)

        refinedMarkerSet = MarkerSet(markerSet.UID, markerSet.lineageStr, markerSet.numGenomes, finalMarkerSet)

        return refinedMarkerSet

    def getBinMarkerSets(self, outDir, markerFile,
                                    numGenomesMarkers,
                                    bootstrap, bNoLineageSpecificRefinement,
                                    bForceDomain, bRequireTaxonomy,
                                    resultsParser, minUnique, maxMulti):
        """Determine marker sets for each bin."""

        self.logger.info('Determining marker sets for each genome bin.')

        # get all bin ids
        binIds = getBinIdsFromOutDir(outDir)

        # get statistics for internal nodes
        uniqueIdToLineageStatistics = self.readNodeMetadata()

        # determine marker set for each bin
        treeFile = os.path.join(outDir, 'storage', 'tree', DefaultValues.PPLACER_TREE_OUT)
        tree = dendropy.Tree.get_from_path(treeFile, schema='newick', rooting="force-rooted", preserve_underscores=True)
        rootNode = tree.find_node(filter_fn=lambda n: n.parent_node == None)

        fout = open(markerFile, 'w')
        fout.write(DefaultValues.LINEAGE_MARKER_FILE_HEADER + '\n')

        numProcessedBins = 0
        statusStr = ''
        for binId in binIds:
            if self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedBins += 1
                sys.stderr.write(' ' * len(statusStr) + '\r')  # clear previous line
                statusStr = '    Finished processing %d of %d (%.2f%%) bins (current: %s).' % (numProcessedBins, len(binIds), float(numProcessedBins) * 100 / len(binIds), binId)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

            node = tree.find_node_with_taxon_label(binId)
            binMarkerSets = BinMarkerSets(binId, BinMarkerSets.TREE_MARKER_SET)
            if node == None:
                # bin is not in tree
                node, markerSet = self.__getMarkerSet(rootNode, tree, uniqueIdToLineageStatistics,
                                                        numGenomesMarkers, bootstrap,
                                                        bForceDomain, bRequireTaxonomy)
                binMarkerSets.addMarkerSet(markerSet)
            else:
                # special case: if node is on the bacterial or archaeal branch descendant from the root,
                # then move down the tree to include the domain-specific marker set
                parentNode = node.parent_node
                while parentNode != None:
                    if parentNode.label:
                        bRoot = (parentNode.parent_node == None)
                        break

                    parentNode = parentNode.parent_node

                if bRoot:
                    # since the root is the first labeled node, we need to descend the
                    # tree to incorporate the domain-specific marker set
                    domainNode = self.__findDomainNode(node)
                    curNode = domainNode.child_nodes()[0]
                else:
                    curNode = node

                # get lineage specific refinement for first node with an id
                if not bNoLineageSpecificRefinement:
                    uniqueId = parentNode.label.split('|')[0]
                    self.__readLineageSpecificGenesToRemove()
                    lineageSpecificRefinement = self.lineageSpecificGenesToRemove[uniqueId]

                # ascend tree to root, recording all marker sets meeting selection criteria
                while curNode.parent_node != None:
                    uniqueHits, multiCopyHits = resultsParser.results[binId].countUniqueHits()
                    tempForceDomain = bForceDomain or (uniqueHits < minUnique) or (multiCopyHits > maxMulti)

                    curNode, markerSet = self.__getMarkerSet(curNode.parent_node, tree, uniqueIdToLineageStatistics,
                                                                numGenomesMarkers, bootstrap,
                                                                tempForceDomain, bRequireTaxonomy)

                    if not bNoLineageSpecificRefinement:
                        markerSet = self.__removeInvalidLineageMarkerGenes(markerSet, lineageSpecificRefinement)

                    binMarkerSets.addMarkerSet(markerSet)

            binMarkerSets.write(fout)

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

        fout.close()

    def readNodeMetadata(self):
        """Read metadata for internal nodes."""

        uniqueIdToLineageStatistics = {}
        metadataFile = os.path.join(DefaultValues.GENOME_TREE_DIR, DefaultValues.GENOME_TREE_METADATA)
        with open(metadataFile) as f:
            f.readline()
            for line in f:
                lineSplit = line.rstrip().split('\t')

                uniqueId = lineSplit[0]

                d = {}
                d['# genomes'] = int(lineSplit[1])
                d['taxonomy'] = lineSplit[2]
                try:
                    d['bootstrap'] = float(lineSplit[3])
                except:
                    d['bootstrap'] = 'NA'
                d['gc mean'] = float(lineSplit[4])
                d['gc std'] = float(lineSplit[5])
                d['genome size mean'] = float(lineSplit[6]) / 1e6
                d['genome size std'] = float(lineSplit[7]) / 1e6
                d['gene count mean'] = float(lineSplit[8])
                d['gene count std'] = float(lineSplit[9])
                d['marker set'] = lineSplit[10].rstrip()

                uniqueIdToLineageStatistics[uniqueId] = d

        return uniqueIdToLineageStatistics

    def readLineageMetadata(self, outDir, binIds):
        """Get metadata for each bin."""

        uniqueIdToLineageStatistics = self.readNodeMetadata()

        # read genome tree
        treeFile = os.path.join(outDir, 'storage', 'tree', DefaultValues.PPLACER_TREE_OUT)
        tree = dendropy.Tree.get_from_path(treeFile, schema='newick', rooting="force-rooted", preserve_underscores=True)

        # find first parent of each bin with a label
        binIdToLineageStatistics = {}
        for binId in binIds:
            node = tree.find_node_with_taxon_label(binId)
            if node == None:
                d = {}
                d['# genomes'] = 'NA'
                d['taxonomy'] = 'unresolved'
                d['gc mean'] = 'NA'
                d['gc std'] = 'NA'
                d['genome size mean'] = 'NA'
                d['genome size std'] = 'NA'
                d['gene count mean'] = 'NA'
                d['gene count std'] = 'NA'
                d['marker set'] = 'NA'
                binIdToLineageStatistics[binId] = d
                continue

            # find first labeled parent node (nodes inserted by pplacer will be unlabeled)
            parentNode = node.parent_node
            uniqueId = None
            while parentNode != None:
                if parentNode.label:
                    uniqueId = parentNode.label.split('|')[0]
                    break

                parentNode = parentNode.parent_node

            if uniqueId:
                binIdToLineageStatistics[binId] = uniqueIdToLineageStatistics[uniqueId]
            else:
                self.logger.error('Failed to find lineage-specific statistics for inserted bin: ' + node.taxon.label)
                sys.exit(1)

        return binIdToLineageStatistics

    def reportBinTaxonomy(self, outDir, resultsParser, bTabTable, outFile, binStats, bLineageStatistics):
        # make sure output and tree directories exist
        checkDirExists(outDir)
        alignOutputDir = os.path.join(outDir, 'storage', 'tree')
        checkDirExists(alignOutputDir)

        # get all bin ids
        binIds = getBinIdsFromOutDir(outDir)

        # get taxonomy for each bin
        binIdToTaxonomy = self.getBinTaxonomy(outDir, binIds)

        # write table
        if not bLineageStatistics:
            self.__printSimpleSummaryTable(binIdToTaxonomy, resultsParser, bTabTable, outFile)
        else:
            # get taxonomy of sister lineage for each bin
            binIdToSisterTaxonomy = self.getBinSisterTaxonomy(outDir, binIds)
            binIdToUID = self.getInsertionBranchId(outDir, binIds)

            binIdToLineageStatistics = self.readLineageMetadata(outDir, binIds)
            self.__printFullTable(binIdToUID, binIdToTaxonomy, binIdToSisterTaxonomy, binIdToLineageStatistics, resultsParser, binStats, bTabTable, outFile)

    def __printSimpleSummaryTable(self, binIdToTaxonomy, resultsParser, bTabTable, outFile):
        # redirect output
        oldStdOut = reassignStdOut(outFile)

        arbitraryBinId = list(binIdToTaxonomy.keys())[0]
        markerCountLabel = '# unique markers (of %d)' % len(resultsParser.models[arbitraryBinId])
        header = ['Bin Id', markerCountLabel, '# multi-copy', 'Taxonomy']

        if bTabTable:
            pTable = None
            print('\t'.join(header))
        else:
            pTable = prettytable.PrettyTable(header)
            pTable.float_format = '.2'
            pTable.align = 'c'
            pTable.align[header[0]] = 'l'
            pTable.align['Taxonomy'] = 'l'
            pTable.hrules = prettytable.FRAME
            pTable.vrules = prettytable.NONE

        for binId in sorted(binIdToTaxonomy.keys()):
            uniqueHits, multiCopyHits = resultsParser.results[binId].countUniqueHits()

            row = [binId, uniqueHits, multiCopyHits, binIdToTaxonomy[binId]]

            if bTabTable:
                print('\t'.join(list(map(str, row))))
            else:
                pTable.add_row(row)

        if not bTabTable:
            print(pTable.get_string(sortby=markerCountLabel, reversesort=True))

        # restore stdout
        restoreStdOut(outFile, oldStdOut)

    def __printFullTable(self, binIdToUID, binIdToTaxonomy, binIdToSisterTaxonomy, binIdToLineageStatistics, resultsParser, binStats, bTabTable, outFile):
        # redirect output
        oldStdOut = reassignStdOut(outFile)

        arbitraryBinId = list(binIdToTaxonomy.keys())[0]
        markerCountLabel = '# unique markers (of %d)' % len(resultsParser.models[arbitraryBinId])
        header = ['Bin Id', markerCountLabel, "# multi-copy"]
        header += ['Insertion branch UID', 'Taxonomy (contained)', 'Taxonomy (sister lineage)']
        header += ['GC', 'Genome size (Mbp)', 'Gene count', 'Coding density', 'Translation table']
        header += ['# descendant genomes', 'Lineage: GC mean', 'Lineage: GC std']
        header += ['Lineage: genome size (Mbp) mean', 'Lineage: genome size (Mbp) std']
        header += ['Lineage: gene count mean', 'Lineage: gene count std']

        if bTabTable:
            pTable = None
            print('\t'.join(header))
        else:
            pTable = prettytable.PrettyTable(header)
            pTable.float_format = '.2'
            pTable.float_format['GC'] = '.1'
            pTable.float_format['Lineage: GC mean'] = '.1'
            pTable.float_format['Lineage: GC std'] = '.1'
            pTable.float_format['Lineage: gene count mean'] = '.0'
            pTable.float_format['Lineage: gene count std'] = '.0'
            pTable.align = 'c'
            pTable.align[header[0]] = 'l'
            pTable.align['Insertion branch UID'] = 'l'
            pTable.align['Taxonomy (contained)'] = 'l'
            pTable.align['Taxonomy (sister lineage)'] = 'l'
            pTable.hrules = prettytable.FRAME
            pTable.vrules = prettytable.NONE

        for binId in sorted(binIdToTaxonomy.keys()):
            uniqueHits, multiCopyHits = resultsParser.results[binId].countUniqueHits()

            truncSisterLineage = binIdToSisterTaxonomy[binId]
            for taxa in binIdToTaxonomy[binId].split(';'):
                truncSisterLineage = truncSisterLineage.replace(taxa + ';', '')

            if len(truncSisterLineage) == 0:
                truncSisterLineage = 'unresolved'
            elif truncSisterLineage[-1] == ';':
                truncSisterLineage = truncSisterLineage[0:-1]

            row = [binId, uniqueHits, multiCopyHits]
            row += [binIdToUID[binId], binIdToTaxonomy[binId], truncSisterLineage]
            row += [binStats[binId]['GC'] * 100]
            row += [float(binStats[binId]['Genome size']) / 1e6]
            row += [binStats[binId]['# predicted genes']]
            row += [binStats[binId]['Coding density']]
            row += [binStats[binId]['Translation table']]
            row += [binIdToLineageStatistics[binId]['# genomes']]
            row += [binIdToLineageStatistics[binId]['gc mean']]
            row += [binIdToLineageStatistics[binId]['gc std']]
            row += [binIdToLineageStatistics[binId]['genome size mean']]
            row += [binIdToLineageStatistics[binId]['genome size std']]
            row += [binIdToLineageStatistics[binId]['gene count mean']]
            row += [binIdToLineageStatistics[binId]['gene count std']]

            if bTabTable:
                print('\t'.join(list(map(str, row))))
            else:
                pTable.add_row(row)

        if not bTabTable:
            print(pTable.get_string(sortby=markerCountLabel, reversesort=True))

        # restore stdout
        restoreStdOut(outFile, oldStdOut)

    def __readDuplicateSeqs(self):
        """Parse file indicating duplicate sequence alignments."""
        duplicateSeqs = {}
        for line in open(os.path.join(DefaultValues.GENOME_TREE_DIR, DefaultValues.GENOME_TREE_DEREP)):
            lineSplit = line.rstrip().split()
            if len(lineSplit) > 1:
                duplicateSeqs[lineSplit[0]] = lineSplit[1:]

        return duplicateSeqs
 class TaxonParser():
    """Parse taxonomic-specific marker sets."""
    def __init__(self):
        self.logger = logging.getLogger('timestamp')

    def readMarkerSets(self):
        taxonMarkerSets = defaultdict(dict)
        for line in open(DefaultValues.TAXON_MARKER_SETS):
            lineSplit = line.split('\t')
            rank = lineSplit[0]
            taxon = lineSplit[1]
            lineage = lineSplit[2]
            numGenomes = int(lineSplit[3])
            markerSet = eval(lineSplit[6].rstrip())
            
            ms = MarkerSet(ranksByLabel[rank], lineage, numGenomes, markerSet)
            ms.removeMarkers(DefaultValues.MARKERS_TO_EXCLUDE)

            taxonMarkerSets[rank][taxon] = ms

        return taxonMarkerSets

    def list(self, rankFilter='ALL'):
        """ List all available marker sets from the specified rank."""

        taxonMarkerSets = self.readMarkerSets()

        header = ['Rank', 'Taxon', '# genomes', '# marker genes', '# marker sets']
        pTable = prettytable.PrettyTable(header)
        pTable.align = 'c'
        pTable.align['Rank'] = 'l'
        pTable.align['Taxon'] = 'l'
        pTable.hrules = prettytable.FRAME
        pTable.vrules = prettytable.NONE

        for rank in taxonomicRanks:
            if rankFilter == 'ALL' or rankFilter == rank:
                for taxon in sorted(taxonMarkerSets[rank]):
                    markerSet = taxonMarkerSets[rank][taxon]

                    numMarkers, numMarkerSets = markerSet.size()
                    pTable.add_row([rank, taxon, markerSet.numGenomes, numMarkers, numMarkerSets])

        print ('')
        print (pTable.get_string())

    def markerSet(self, rank, taxon, markerFile):
        """Obtain specified taxonomic-specific marker set."""

        taxonMarkerSets = self.readMarkerSets()

        if rank not in taxonMarkerSets:
            self.logger.error('Unrecognized taxonomic rank: ' + rank)
            return False
        elif taxon not in taxonMarkerSets[rank]:
            self.logger.error('Unrecognized taxon: %s (in rank %s): ' % (taxon, rank))
            return False

        markerSet = taxonMarkerSets[rank][taxon]

        taxonomy = markerSet.lineageStr.split(';')[::-1]
        binMarkerSets = BinMarkerSets(taxon, BinMarkerSets.TAXONOMIC_MARKER_SET)
        for i, taxon in enumerate(taxonomy):
            if rank != 'life':
                rank = ranksByLevel[len(taxonomy) - i - 1]

            if rank == 'species':
                taxon = taxonomy[1] + ' ' + taxonomy[0]

            markerSet = taxonMarkerSets[rank][taxon]
            numMarkers, numMarkerSets = markerSet.size()
            self.logger.info('Marker set for %s contains %d marker genes arranged in %d sets.' % (taxon, numMarkers, numMarkerSets))
            self.logger.info('Marker set inferred from %d reference genomes.' % markerSet.numGenomes)

            markerSet.lineageStr = taxon
            binMarkerSets.addMarkerSet(markerSet)

        fout = open(markerFile, 'w')
        fout.write(DefaultValues.TAXON_MARKER_FILE_HEADER + '\n')
        binMarkerSets.write(fout)
        fout.close()

        return True
 class AminoAcidIdentity():
    """Calculate AAI between sequences aligned to an HMM."""
    def __init__(self):
        self.logger = logging.getLogger('timestamp')
        self.aaiRawScores = defaultdict(dict)
        self.aaiHetero = defaultdict(dict)
        self.aaiMeanBinHetero = {}

    def run(self, aaiStrainThreshold, outDir, alignmentOutputFile):
        """Calculate AAI between input alignments."""

        self.logger.info('Calculating AAI between multi-copy marker genes.')

        if alignmentOutputFile:
            fout = open(alignmentOutputFile, 'w')

        # calculate AAI for duplicate marker genes
        binIds = getBinIdsFromOutDir(outDir)
        aaiOutputDir = os.path.join(outDir, 'storage', 'aai_qa')
        for binId in binIds:
            binPath = os.path.join(aaiOutputDir, binId)
            if not os.path.exists(binPath):
                continue

            for f in os.listdir(binPath):
                if not f.endswith('.masked.faa'):
                    continue

                markerId = f[0:f.find('.')]

                seqs = readFasta(os.path.join(binPath, f))

                # calculate AAI between all pairs of seqs
                for i in range(0, len(seqs)):
                    seqIdI = list(seqs.keys())[i]
                    binIdI = seqIdI[0:seqIdI.find(DefaultValues.SEQ_CONCAT_CHAR)]

                    seqI = seqs[seqIdI]

                    for j in range(i + 1, len(seqs)):
                        seqIdJ = list(seqs.keys())[j]
                        binIdJ = seqIdJ[0:seqIdJ.find(DefaultValues.SEQ_CONCAT_CHAR)]

                        seqJ = seqs[seqIdJ]

                        if binIdI == binIdJ:
                            aai = self.aai(seqI, seqJ)

                            if alignmentOutputFile:
                                fout.write(binId + ',' + markerId + '\n')
                                fout.write(seqIdI + '\t' + seqI + '\n')
                                fout.write(seqIdJ + '\t' + seqJ + '\n')
                                fout.write('AAI: %.3f\n' % aai)
                                fout.write('\n')

                            if binIdI not in self.aaiRawScores:
                                self.aaiRawScores[binIdI] = defaultdict(list)
                            self.aaiRawScores[binIdI][markerId].append(aai)
                        else:
                            # something is wrong as the bin Ids should always be the same
                            self.logger.error('Bin ids do not match.')
                            sys.exit(1)

        if alignmentOutputFile:
            fout.close()

        # calculate strain heterogeneity for each marker gene in each bin
        self.aaiHetero, self.aaiMeanBinHetero = self.strainHetero(self.aaiRawScores, aaiStrainThreshold)

    def strainHetero(self, aaiScores, aaiStrainThreshold):
        """Calculate strain heterogeneity."""
        aaiHetero = defaultdict(dict)
        aaiMeanBinHetero = {}

        for binId, markerIds in aaiScores.items():
            strainCount = 0
            multiCopyPairs = 0

            aaiHetero[binId] = {}

            for markerId, aaiScores in markerIds.items():
                localStrainCount = 0
                for aaiScore in aaiScores:
                    multiCopyPairs += 1
                    if aaiScore > aaiStrainThreshold:
                        strainCount += 1
                        localStrainCount += 1

                strainHetero = float(localStrainCount) / len(aaiScores)
                aaiHetero[binId][markerId] = strainHetero

            aaiMeanBinHetero[binId] = 100 * float(strainCount) / multiCopyPairs

        return aaiHetero, aaiMeanBinHetero

    def aai(self, seq1, seq2):
        """Calculate amino acid identity between sequences."""
        assert len(seq1) == len(seq2)

        # calculation of AAI should ignore missing data at
        # the start of end of each sequence
        startIndex = 0
        for i in range(0, len(seq1)):
            if seq1[i] == '-' or seq2[i] == '-':
                startIndex = i + 1
            else:
                break

        endIndex = len(seq1)
        for i in range(len(seq1) - 1, 0, -1):
            if seq1[i] == '-' or seq2[i] == '-':
                endIndex = i
            else:
                break

        mismatches = 0
        seqLen = 0
        for i in range(startIndex, endIndex):
            if seq1[i] != seq2[i]:
                mismatches += 1
                seqLen += 1
            elif seq1[i] == '-' and seq2[i] == '-':
                pass
            else:
                seqLen += 1

        if seqLen == 0:
            return 0.0

        return 1.0 - (float(mismatches) / seqLen)
 class BinStatistics():
    """Statistics for genomes.

    This class calculates statistics for genomes comprised
    of one or more scaffolds. Mean statistics are weighted by
    scaffold length. The following statistics are calculated:
     - bin assignment
     - mean GC
     - mean and median scaffold length
     - N50
     - mean coverage
     - mean tetranucleotide signature

    """

    def __init__(self, threads=1):
        """Initialization.

        Parameters
        ----------
        cpus : int
            Number of cpus to use.
        """

        self.logger = logging.getLogger('timestamp')
        self.totalThreads = threads

    def calculate(self, binFiles, outDir, binStatsFile):
        """Calculate statistics for each putative genome bin."""

        # process each bin
        self.logger.info("Calculating genome statistics for %d bins with %d threads:" % (len(binFiles), self.totalThreads))

        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        for binFile in binFiles:
            workerQueue.put(binFile)

        for _ in range(self.totalThreads):
            workerQueue.put(None)

        try:
            calcProc = [mp.Process(target=self.__processBin, args=(outDir, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__reportProgress, args=(outDir, binStatsFile, len(binFiles), writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put((None, None))
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

    def __processBin(self, outDir, queueIn, queueOut):
        """Thread safe bin processing."""
        while True:
            binFile = queueIn.get(block=True, timeout=None)
            if binFile == None:
                break

            binStats = {}

            binId = binIdFromFilename(binFile)
            binDir = os.path.join(outDir, 'bins', binId)
            makeSurePathExists(binDir)

            # read scaffolds
            scaffolds = readFasta(binFile)

            # calculate GC statistics
            GC, stdGC = self.calculateGC(scaffolds)
            binStats['GC'] = GC
            binStats['GC std'] = stdGC

            # calculate statistics related to contigs and scaffolds
            maxScaffoldLen, maxContigLen, genomeSize, scaffold_N50, contig_N50, scaffoldAvgLen, contigAvgLen, numContigs, numAmbiguousBases = self.calculateSeqStats(scaffolds)
            binStats['Genome size'] = genomeSize
            binStats['# ambiguous bases'] = numAmbiguousBases
            binStats['# scaffolds'] = len(scaffolds)
            binStats['# contigs'] = numContigs
            binStats['Longest scaffold'] = maxScaffoldLen
            binStats['Longest contig'] = maxContigLen
            binStats['N50 (scaffolds)'] = scaffold_N50
            binStats['N50 (contigs)'] = contig_N50
            binStats['Mean scaffold length'] = scaffoldAvgLen
            binStats['Mean contig length'] = contigAvgLen

            # calculate coding density statistics
            codingDensity, translationTable, numORFs = self.calculateCodingDensity(binDir, scaffolds, genomeSize)
            binStats['Coding density'] = codingDensity
            binStats['Translation table'] = translationTable
            binStats['# predicted genes'] = numORFs

            queueOut.put((binId, binStats))

    def __reportProgress(self, outDir, binStatsFile, numBins, queueIn):
        """Report number of processed bins and write statistics to file."""

        storagePath = os.path.join(outDir, 'storage')
        fout = open(os.path.join(storagePath, binStatsFile), 'w')

        numProcessedBins = 0
        if self.logger.getEffectiveLevel() <= logging.INFO:
            statusStr = '    Finished processing %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
            sys.stderr.write('%s\r' % statusStr)
            sys.stderr.flush()

        while True:
            binId, curBinStats = queueIn.get(block=True, timeout=None)
            if binId == None:
                break

            fout.write(binId + '\t' + str(curBinStats) + '\n')
            if self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedBins += 1
                statusStr = '    Finished processing %d of %d (%.2f%%) bins.' % (numProcessedBins, numBins, float(numProcessedBins) * 100 / numBins)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

        fout.close()

    def calculateGC(self, seqs, seqStats=None):
        """Calculate fraction of nucleotides that are G or C."""
        totalGC = 0
        totalAT = 0
        gcPerSeq = []
        for seqId, seq in seqs.items():
            a, c, g, t = baseCount(seq)

            gc = g + c
            at = a + t

            totalGC += gc
            totalAT += at

            if (gc + at) > 0:
                gcContent = float(gc) / (gc + at)
            else:
                gcContent = 0.0

            if seqStats:
                seqStats[seqId]['GC'] = gcContent

            if len(seq) > DefaultValues.MIN_SEQ_LEN_GC_STD:
                gcPerSeq.append(gcContent)

        if (totalGC + totalAT) > 0:
            GC = float(totalGC) / (totalGC + totalAT)
        else:
            GC = 0.0

        varGC = 0
        if len(gcPerSeq) > 1:
            varGC = mean(list(map(lambda x: (x - GC) ** 2, gcPerSeq)))

        return GC, math.sqrt(varGC)

    def calculateSeqStats(self, scaffolds, seqStats=None):
        """Calculate scaffold length statistics (min length, max length, total length, N50, # contigs)."""
        scaffoldLens = []
        contigLens = []
        numAmbiguousBases = 0
        for scaffoldId, scaffold in scaffolds.items():
            scaffoldLen = len(scaffold)
            scaffoldLens.append(scaffoldLen)

            splitScaffold = scaffold.split(DefaultValues.CONTIG_BREAK)
            lenContigsInScaffold = []
            for contig in splitScaffold:
                contigLen = len(contig.replace('N', ''))
                if contigLen > 0:
                    lenContigsInScaffold.append(contigLen)

            contigLens += lenContigsInScaffold

            if seqStats:
                seqStats[scaffoldId]['Length'] = scaffoldLen
                seqStats[scaffoldId]['Total contig length'] = sum(lenContigsInScaffold)
                seqStats[scaffoldId]['# contigs'] = len(lenContigsInScaffold)

            numAmbiguousBases += scaffold.count('N') + scaffold.count('n')

        scaffold_N50 = calculateN50(scaffoldLens)
        contig_N50 = calculateN50(contigLens)

        return max(scaffoldLens), max(contigLens), sum(scaffoldLens), scaffold_N50, contig_N50, mean(scaffoldLens), mean(contigLens), len(contigLens), numAmbiguousBases

    def calculateCodingDensity(self, outDir, scaffolds, genomeSize):
        """Calculate coding density of putative genome bin."""
        gffFile = os.path.join(outDir, DefaultValues.PRODIGAL_GFF)
        if os.path.exists(gffFile):
            prodigalParserGFF = ProdigalGeneFeatureParser(gffFile)

            aaFile = os.path.join(outDir, DefaultValues.PRODIGAL_AA)  # use AA file as nucleotide file is optional
            aaGenes = readFasta(aaFile)

            codingBasePairs = 0  # self.__calculateCodingBases(aaGenes)
            for scaffold_id in scaffolds.keys():
                codingBasePairs += prodigalParserGFF.codingBases(scaffold_id)

            return float(codingBasePairs) / genomeSize, prodigalParserGFF.translationTable, len(aaGenes)
        else:
            # there is no gene feature file (perhaps the user specified pre-calculated genes)
            # so calculating the coding density is not possible
            return -1, -1, -1

    def __calculateCodingBases(self, aaGenes):
        """Calculate number of coding bases in a set of genes."""
        codingBasePairs = 0
        for _geneId, gene in aaGenes.items():
            codingBasePairs += len(gene) * 3

        return codingBasePairs

    def sequenceStats(self, outDir, binFile):
        """Calculate statistics for all sequences within a bin."""

        # read scaffolds
        seqs = readFasta(binFile)

        seqStats = {}
        for seqId in seqs:
            seqStats[seqId] = {}

        self.calculateGC(seqs, seqStats)
        self.calculateSeqStats(seqs, seqStats)

        binId = binIdFromFilename(binFile)
        aaFile = os.path.join(outDir, 'bins', binId, DefaultValues.PRODIGAL_AA)
        if os.path.exists(aaFile):
            aaGenes = readFasta(aaFile)
            for geneId, gene in aaGenes.items():
                seqId = geneId[0:geneId.rfind('_')]
                seqStats[seqId]['# ORFs'] = seqStats[seqId].get('# ORFs', 0) + 1
                seqStats[seqId]['Coding bases'] = seqStats[seqId].get('Coding bases', 0) + len(gene) * 3
        else:
            # missing amino acid file likely indicates users used a pre-called gene file, so
            # just set some defaults
            seqStats[seqId]['# ORFs'] = seqStats[seqId].get('# ORFs', 0) + 1
            seqStats[seqId]['Coding bases'] = seqStats[seqId].get('Coding bases', 0) + len(gene) * 3

        return seqStats class Coverage():
    """Calculate coverage of all sequences."""
    def __init__(self, threads):
        self.logger = logging.getLogger('timestamp')

        self.totalThreads = threads

    def run(self, binFiles, bamFiles, outFile, bAllReads, minAlignPer, maxEditDistPer, minQC):
        """Calculate coverage of sequences for each BAM file."""

        # determine bin assignment of each sequence
        self.logger.info('Determining bin assignment of each sequence.')

        seqIdToBinId = {}
        seqIdToSeqLen = {}
        for binFile in binFiles:
            binId = binIdFromFilename(binFile)

            seqs = readFasta(binFile)
            for seqId, seq in seqs.items():
                seqIdToBinId[seqId] = binId
                seqIdToSeqLen[seqId] = len(seq)

        # process each fasta file
        self.logger.info("Processing %d file(s) with %d threads.\n" % (len(bamFiles), self.totalThreads))

        # make sure all BAM files are sorted
        self.numFiles = len(bamFiles)
        for bamFile in bamFiles:
            if not os.path.exists(bamFile + '.bai'):
                self.logger.error('BAM file is either unsorted or not indexed: ' + bamFile + '\n')
                sys.exit(1)

        # calculate coverage of each BAM file
        coverageInfo = {}
        numFilesStarted = 0
        for bamFile in bamFiles:
            numFilesStarted += 1
            self.logger.info('Processing %s (%d of %d):' % (ntpath.basename(bamFile), numFilesStarted, len(bamFiles)))

            coverageInfo[bamFile] = mp.Manager().dict()
            coverageInfo[bamFile] = self.__processBam(bamFile, bAllReads, minAlignPer, maxEditDistPer, minQC, coverageInfo[bamFile])

        # redirect output
        self.logger.info('Writing coverage information to file.')
        oldStdOut = reassignStdOut(outFile)

        header = 'Sequence Id\tBin Id\tSequence length (bp)'
        for bamFile in bamFiles:
            header += '\tBam Id\tCoverage\tMapped reads'

        print(header)

        # get length of all seqs
        for bamFile, seqIds in coverageInfo.items():
            for seqId in seqIds.keys():
                seqIdToSeqLen[seqId] = seqIds[seqId].seqLen

        # write coverage stats for all scaffolds to file
        for seqId, seqLen in seqIdToSeqLen.items():
            rowStr = seqId + '\t' + seqIdToBinId.get(seqId, DefaultValues.UNBINNED) + '\t' + str(seqLen)
            for bamFile in bamFiles:
                bamId = binIdFromFilename(bamFile)

                if seqId in coverageInfo[bamFile]:
                    rowStr += '\t%s\t%f\t%d' % (bamId, coverageInfo[bamFile][seqId].coverage, coverageInfo[bamFile][seqId].mappedReads)
                else:
                    rowStr += '\t%s\t%f\t%d' % (bamId, 0, 0)

            print(rowStr)

        # restore stdout
        restoreStdOut(outFile, oldStdOut)

    def __processBam(self, bamFile, bAllReads, minAlignPer, maxEditDistPer, minQC, coverageInfo):
        """Calculate coverage of sequences in BAM file."""

        # determine coverage for each reference sequence
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        bamfile = pysam.Samfile(bamFile, 'rb')
        refSeqIds = bamfile.references
        refSeqLens = bamfile.lengths

        # populate each thread with reference sequence to process
        # Note: reference sequences are sorted by number of mapped reads
        # so it is important to distribute reads in a sensible way to each
        # of the threads
        refSeqLists = [[] for _ in range(self.totalThreads)]
        refLenLists = [[] for _ in range(self.totalThreads)]

        threadIndex = 0
        incDir = 1
        for refSeqId, refLen in zip(refSeqIds, refSeqLens):
            refSeqLists[threadIndex].append(refSeqId)
            refLenLists[threadIndex].append(refLen)

            threadIndex += incDir
            if threadIndex == self.totalThreads:
                threadIndex = self.totalThreads - 1
                incDir = -1
            elif threadIndex == -1:
                threadIndex = 0
                incDir = 1

        for i in range(self.totalThreads):
            workerQueue.put((refSeqLists[i], refLenLists[i]))

        for _ in range(self.totalThreads):
            workerQueue.put((None, None))

        try:
            workerProc = [mp.Process(target=self.__workerThread, args=(bamFile, bAllReads, minAlignPer, maxEditDistPer, minQC, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__writerThread, args=(coverageInfo, len(refSeqIds), writerQueue))

            writeProc.start()

            for p in workerProc:
                p.start()

            for p in workerProc:
                p.join()

            writerQueue.put((None, None, None, None, None, None, None, None, None, None, None))
            writeProc.join()
        except:
            # make sure all processes are terminated
            print (traceback.format_exc())
            for p in workerProc:
                p.terminate()

            writeProc.terminate()

        return coverageInfo

    def __workerThread(self, bamFile, bAllReads, minAlignPer, maxEditDistPer, minQC, queueIn, queueOut):
        """Process each data item in parallel."""
        while True:
            seqIds, seqLens = queueIn.get(block=True, timeout=None)
            if seqIds == None:
                break

            bamfile = pysam.Samfile(bamFile, 'rb')

            for seqId, seqLen in zip(seqIds, seqLens):
                numReads = 0
                numMappedReads = 0
                numDuplicates = 0
                numSecondary = 0
                numFailedQC = 0
                numFailedAlignLen = 0
                numFailedEditDist = 0
                numFailedProperPair = 0
                coverage = 0

                for read in bamfile.fetch(seqId, 0, seqLen):
                    numReads += 1

                    if read.is_unmapped:
                        pass
                    elif read.is_duplicate:
                        numDuplicates += 1
                    elif read.is_secondary or read.is_supplementary:
                        numSecondary += 1
                    elif read.is_qcfail or read.mapping_quality < minQC:
                        numFailedQC += 1
                    elif read.query_alignment_length < minAlignPer * read.query_length:
                        numFailedAlignLen += 1
                    elif read.get_tag('NM') > maxEditDistPer * read.query_length:
                        numFailedEditDist += 1
                    elif not bAllReads and not read.is_proper_pair:
                        numFailedProperPair += 1
                    else:
                        numMappedReads += 1

                        # Note: the alignment length (query_alignment_length) is used instead of the
                        # read length (query_length) as this bring the calculated coverage
                        # in line with 'samtools depth' (at least when the min
                        # alignment length and edit distance thresholds are zero).
                        coverage += read.query_alignment_length

                coverage = float(coverage) / seqLen

                queueOut.put((seqId, seqLen, coverage, numReads,
                                numDuplicates, numSecondary, numFailedQC,
                                numFailedAlignLen, numFailedEditDist,
                                numFailedProperPair, numMappedReads))

            bamfile.close()

    def __writerThread(self, coverageInfo, numRefSeqs, writerQueue):
        """Store or write results of worker threads in a single thread."""
        totalReads = 0
        totalDuplicates = 0
        totalSecondary = 0
        totalFailedQC = 0
        totalFailedAlignLen = 0
        totalFailedEditDist = 0
        totalFailedProperPair = 0
        totalMappedReads = 0

        processedRefSeqs = 0
        while True:
            seqId, seqLen, coverage, numReads, numDuplicates, numSecondary, numFailedQC, numFailedAlignLen, numFailedEditDist, numFailedProperPair, numMappedReads = writerQueue.get(block=True, timeout=None)
            if seqId == None:
                break

            if self.logger.getEffectiveLevel() <= logging.INFO:
                processedRefSeqs += 1
                statusStr = '    Finished processing %d of %d (%.2f%%) reference sequences.' % (processedRefSeqs, numRefSeqs, float(processedRefSeqs) * 100 / numRefSeqs)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

                totalReads += numReads
                totalDuplicates += numDuplicates
                totalSecondary += numSecondary
                totalFailedQC += numFailedQC
                totalFailedAlignLen += numFailedAlignLen
                totalFailedEditDist += numFailedEditDist
                totalFailedProperPair += numFailedProperPair
                totalMappedReads += numMappedReads

            coverageInfo[seqId] = CoverageStruct(seqLen=seqLen, mappedReads=numMappedReads, coverage=coverage)

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

            print ('')
            print ('    # total reads: %d' % totalReads)
            print ('      # properly mapped reads: %d (%.1f%%)' % (totalMappedReads, float(totalMappedReads) * 100 / totalReads))
            print ('      # duplicate reads: %d (%.1f%%)' % (totalDuplicates, float(totalDuplicates) * 100 / totalReads))
            print ('      # secondary reads: %d (%.1f%%)' % (totalSecondary, float(totalSecondary) * 100 / totalReads))
            print ('      # reads failing QC: %d (%.1f%%)' % (totalFailedQC, float(totalFailedQC) * 100 / totalReads))
            print ('      # reads failing alignment length: %d (%.1f%%)' % (totalFailedAlignLen, float(totalFailedAlignLen) * 100 / totalReads))
            print ('      # reads failing edit distance: %d (%.1f%%)' % (totalFailedEditDist, float(totalFailedEditDist) * 100 / totalReads))
            print ('      # reads not properly paired: %d (%.1f%%)' % (totalFailedProperPair, float(totalFailedProperPair) * 100 / totalReads))
            print ('')

    def parseCoverage(self, coverageFile):
        """Read coverage information from file."""
        coverageStats = {}
        bHeader = True
        for line in open(coverageFile):
            if bHeader:
                bHeader = False
                continue

            lineSplit = line.split('\t')
            seqId = lineSplit[0]
            binId = lineSplit[1]

            if binId not in coverageStats:
                coverageStats[binId] = {}

            if seqId not in coverageStats[binId]:
                coverageStats[binId][seqId] = {}

            for i in range(3, len(lineSplit), 3):
                bamId = lineSplit[i]
                coverage = float(lineSplit[i + 1])
                coverageStats[binId][seqId][bamId] = coverage

        return coverageStats

    def binProfiles(self, coverageFile):
        """Read coverage information for each bin."""
        binCoverages = defaultdict(lambda: defaultdict(list))
        binStats = defaultdict(dict)

        bHeader = True
        for line in open(coverageFile):
            if bHeader:
                bHeader = False
                continue

            lineSplit = line.split('\t')
            binId = lineSplit[1]
            seqLen = int(lineSplit[2])

            # calculate mean coverage (weighted by scaffold length)
            # for each bin under each BAM file
            for i in range(3, len(lineSplit), 3):
                bamId = lineSplit[i]
                coverage = float(lineSplit[i + 1])
                binCoverages[binId][bamId].append(coverage)

                if bamId not in binStats[binId]:
                    binStats[binId][bamId] = [0, 0]

                binLength = binStats[binId][bamId][0] + seqLen
                weight = float(seqLen) / binLength
                meanBinCoverage = coverage * weight + binStats[binId][bamId][1] * (1 - weight)

                binStats[binId][bamId] = [binLength, meanBinCoverage]

        profiles = defaultdict(dict)
        for binId in binStats:
            for bamId, stats in binStats[binId].items():
                binLength, meanBinCoverage = stats
                coverages = binCoverages[binId][bamId]

                varCoverage = 0
                if len(coverages) > 1:
                    varCoverage = mean(list(map(lambda x: (x - meanBinCoverage) ** 2, coverages)))

                profiles[binId][bamId] = [meanBinCoverage, sqrt(varCoverage)]

        return profiles class CoverageWindows():
    """Calculate coverage of all sequences."""
    def __init__(self, threads):
        self.logger = logging.getLogger('timestamp')

        self.totalThreads = threads

    def run(self, binFiles, bamFile, bAllReads, minAlignPer, maxEditDistPer, windowSize):
        """Calculate coverage of full sequences and windows."""

        # make sure BAM file is sorted
        if not os.path.exists(bamFile + '.bai'):
            self.logger.error('BAM file is not sorted: ' + bamFile + '\n')
            sys.exit(1)

        # calculate coverage of each BAM file
        self.logger.info('Calculating coverage of windows.')
        coverageInfo = mp.Manager().dict()
        coverageInfo = self.__processBam(bamFile, bAllReads, minAlignPer, maxEditDistPer, windowSize, coverageInfo)

        return coverageInfo

    def __processBam(self, bamFile, bAllReads, minAlignPer, maxEditDistPer, windowSize, coverageInfo):
        """Calculate coverage of sequences in BAM file."""

        # determine coverage for each reference sequence
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        bamfile = pysam.Samfile(bamFile, 'rb')
        refSeqIds = bamfile.references
        refSeqLens = bamfile.lengths

        # populate each thread with reference sequence to process
        # Note: reference sequences are sorted by number of mapped reads
        # so it is important to distribute reads in a sensible way to each
        # of the threads
        refSeqLists = [[] for _ in range(self.totalThreads)]
        refLenLists = [[] for _ in range(self.totalThreads)]

        threadIndex = 0
        incDir = 1
        for refSeqId, refLen in zip(refSeqIds, refSeqLens):
            refSeqLists[threadIndex].append(refSeqId)
            refLenLists[threadIndex].append(refLen)

            threadIndex += incDir
            if threadIndex == self.totalThreads:
                threadIndex = self.totalThreads - 1
                incDir = -1
            elif threadIndex == -1:
                threadIndex = 0
                incDir = 1

        for i in range(self.totalThreads):
            workerQueue.put((refSeqLists[i], refLenLists[i]))

        for _ in range(self.totalThreads):
            workerQueue.put((None, None))

        try:
            workerProc = [mp.Process(target=self.__workerThread, args=(bamFile, bAllReads, minAlignPer, maxEditDistPer, windowSize, workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__writerThread, args=(coverageInfo, len(refSeqIds), writerQueue))

            writeProc.start()

            for p in workerProc:
                p.start()

            for p in workerProc:
                p.join()

            writerQueue.put((None, None, None, None, None, None, None, None, None, None, None, None))
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in workerProc:
                p.terminate()

            writeProc.terminate()

        return coverageInfo

    def __workerThread(self, bamFile, bAllReads, minAlignPer, maxEditDistPer, windowSize, queueIn, queueOut):
        """Process each data item in parallel."""
        while True:
            seqIds, seqLens = queueIn.get(block=True, timeout=None)
            if seqIds == None:
                break

            bamfile = pysam.Samfile(bamFile, 'rb')

            for seqId, seqLen in zip(seqIds, seqLens):
                readLoader = ReadLoader(seqLen, bAllReads, minAlignPer, maxEditDistPer)
                bamfile.fetch(seqId, 0, seqLen, callback=readLoader)

                start = 0
                end = windowSize
                windowCoverages = []
                while(end < seqLen):
                    windowCoverages.append(sum(readLoader.coverage[start:end]) / windowSize)

                    start = end
                    try:
                        end += windowSize
                    except:
                        print ('*****************')
                        print (end)
                        print (windowSize)
                        print ('******************')

                coverage = float(sum(readLoader.coverage)) / seqLen

                queueOut.put((seqId, seqLen, coverage, windowCoverages, readLoader.numReads,
                                readLoader.numDuplicates, readLoader.numSecondary, readLoader.numFailedQC,
                                readLoader.numFailedAlignLen, readLoader.numFailedEditDist,
                                readLoader.numFailedProperPair, readLoader.numMappedReads))

            bamfile.close()

    def __writerThread(self, coverageInfo, numRefSeqs, writerQueue):
        """Store or write results of worker threads in a single thread."""
        totalReads = 0
        totalDuplicates = 0
        totalSecondary = 0
        totalFailedQC = 0
        totalFailedAlignLen = 0
        totalFailedEditDist = 0
        totalFailedProperPair = 0
        totalMappedReads = 0

        processedRefSeqs = 0
        while True:
            seqId, seqLen, coverage, windowCoverages, numReads, numDuplicates, numSecondary, numFailedQC, numFailedAlignLen, numFailedEditDist, numFailedProperPair, numMappedReads = writerQueue.get(block=True, timeout=None)
            if seqId == None:
                break

            if self.logger.getEffectiveLevel() <= logging.INFO:
                processedRefSeqs += 1
                statusStr = '    Finished processing %d of %d (%.2f%%) reference sequences.' % (processedRefSeqs, numRefSeqs, float(processedRefSeqs) * 100 / numRefSeqs)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

                totalReads += numReads
                totalDuplicates += numDuplicates
                totalSecondary += numSecondary
                totalFailedQC += numFailedQC
                totalFailedAlignLen += numFailedAlignLen
                totalFailedEditDist += numFailedEditDist
                totalFailedProperPair += numFailedProperPair
                totalMappedReads += numMappedReads

            coverageInfo[seqId] = [coverage, windowCoverages]

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

            print ('')
            print ('    # total reads: %d' % totalReads)
            print ('      # properly mapped reads: %d (%.1f%%)' % (totalMappedReads, float(totalMappedReads) * 100 / totalReads))
            print ('      # duplicate reads: %d (%.1f%%)' % (totalDuplicates, float(totalDuplicates) * 100 / totalReads))
            print ('      # secondary reads: %d (%.1f%%)' % (totalSecondary, float(totalSecondary) * 100 / totalReads))
            print ('      # reads failing QC: %d (%.1f%%)' % (totalFailedQC, float(totalFailedQC) * 100 / totalReads))
            print ('      # reads failing alignment length: %d (%.1f%%)' % (totalFailedAlignLen, float(totalFailedAlignLen) * 100 / totalReads))
            print ('      # reads failing edit distance: %d (%.1f%%)' % (totalFailedEditDist, float(totalFailedEditDist) * 100 / totalReads))
            print ('      # reads not properly paired: %d (%.1f%%)' % (totalFailedProperPair, float(totalFailedProperPair) * 100 / totalReads))
            print ('') class GenomicSignatures(object):
    def __init__(self, K, threads):
        self.logger = logging.getLogger('timestamp')

        self.K = K
        self.compl = str.maketrans('ACGT', 'TGCA')
        self.kmerCols, self.kmerToCanonicalIndex = self.__makeKmerColNames()

        self.totalThreads = threads

    def __makeKmerColNames(self):
        """Work out unique kmers."""

        # determine all mers of a given length
        baseWords = ("A", "C", "G", "T")
        mers = ["A", "C", "G", "T"]
        for _ in range(1, self.K):
            workingList = []
            for mer in mers:
                for char in baseWords:
                    workingList.append(mer + char)
            mers = workingList

        # pare down kmers based on lexicographical ordering
        retList = []
        for mer in mers:
            kmer = self.__lexicographicallyLowest(mer)
            if kmer not in retList:
                retList.append(kmer)

        sorted(retList)

        # create mapping from kmers to their canonical order position
        kmerToCanonicalIndex = {}
        for index, kmer in enumerate(retList):
            kmerToCanonicalIndex[kmer] = index
            kmerToCanonicalIndex[self.__revComp(kmer)] = index

        return retList, kmerToCanonicalIndex

    def __lexicographicallyLowest(self, seq):
        """Return the lexicographically lowest form of this sequence."""
        rseq = self.__revComp(seq)
        if(seq < rseq):
            return seq
        return rseq

    def __revComp(self, seq):
        """Return the reverse complement of a sequence."""
        # build a dictionary to know what letter to switch to
        return seq.translate(self.compl)[::-1]

    def __calculateResults(self, queueIn, queueOut):
        """Calculate genomic signature of sequences in parallel."""
        while True:
            seqId, seq = queueIn.get(block=True, timeout=None)
            if seqId == None:
                break

            sig = self.seqSignature(seq)

            queueOut.put((seqId, sig))

    def __storeResults(self, seqFile, outputFile, totalSeqs, writerQueue):
        """Store genomic signatures to file."""

        # write header
        fout = open(outputFile, 'w')
        fout.write('Sequence Id')
        for kmer in self.canonicalKmerOrder():
            fout.write('\t' + kmer)
        fout.write('\n')

        numProcessedSeq = 0
        while True:
            seqId, sig = writerQueue.get(block=True, timeout=None)
            if seqId == None:
                break

            if self.logger.getEffectiveLevel() <= logging.INFO:
                numProcessedSeq += 1
                statusStr = '    Finished processing %d of %d (%.2f%%) sequences.' % (numProcessedSeq, totalSeqs, float(numProcessedSeq) * 100 / totalSeqs)
                sys.stderr.write('%s\r' % statusStr)
                sys.stderr.flush()

            fout.write(seqId)
            fout.write('\t' + '\t'.join(list(map(str, sig))))
            fout.write('\n')

        if self.logger.getEffectiveLevel() <= logging.INFO:
            sys.stderr.write('\n')

        fout.close()

    def canonicalKmerOrder(self):
        return self.kmerCols

    def seqSignature(self, seq):
        sig = [0] * len(self.kmerCols)

        tmp_seq = seq.upper()

        numMers = len(tmp_seq) - self.K + 1
        for i in range(0, numMers):
            try:
                kmerIndex = self.kmerToCanonicalIndex[tmp_seq[i:i + self.K]]
                sig[kmerIndex] += 1  # Note: a numpy array would be slow here due to this single element increment
            except KeyError:
                # unknown kmer (e.g., contains a N)
                pass

        # normalize
        sig = np.array(sig, dtype=float)
        sig /= np.sum(sig)

        return sig

    def calculate(self, seqFile, outputFile):
        """Calculate genomic signature of each sequence."""

        self.logger.info('Determining tetranucleotide signature of each sequence.')

        # process each sequence in parallel
        workerQueue = mp.Queue()
        writerQueue = mp.Queue()

        seqs = readFasta(seqFile)

        for seqId, seq in seqs.items():
            workerQueue.put((seqId, seq))

        for _ in range(self.totalThreads):
            workerQueue.put((None, None))

        try:
            calcProc = [mp.Process(target=self.__calculateResults, args=(workerQueue, writerQueue)) for _ in range(self.totalThreads)]
            writeProc = mp.Process(target=self.__storeResults, args=(seqFile, outputFile, len(seqs), writerQueue))

            writeProc.start()

            for p in calcProc:
                p.start()

            for p in calcProc:
                p.join()

            writerQueue.put((None, None))
            writeProc.join()
        except:
            # make sure all processes are terminated
            for p in calcProc:
                p.terminate()

            writeProc.terminate()

    def distance(self, sig1, sig2):
        return np.sum(np.abs(sig1 - sig2))

    def read(self, tetraProfileFile):
        sig = {}
        with open(tetraProfileFile) as f:
            next(f)
            for line in f:
                lineSplit = line.split('\t')
                sig[lineSplit[0]] = np.array([float(x) for x in lineSplit[1:]])

        return sig class Unbinned():
    def __init__(self):
        self.logger = logging.getLogger('timestamp')

    def run(self, binFiles, seqFile, outSeqFile, outStatsFile, minSeqLen):
        checkFileExists(seqFile)

        # get list of sequences in bins
        self.logger.info('Reading binned sequences.')

        binnedSeqs = {}
        totalBinnedBases = 0
        for binFile in binFiles:
            seqs = readFasta(binFile)
            binnedSeqs.update(seqs)
            for seq in seqs.values():
                totalBinnedBases += len(seq)

        self.logger.info('  Read %d (%.2f Mbp) binned sequences.' % (len(binnedSeqs), float(totalBinnedBases) / 1e6))

        # get list of all sequences
        self.logger.info('Reading all sequences.')
        allSeqs = readFasta(seqFile)
        totalBases = 0
        for seq in allSeqs.values():
            totalBases += len(seq)
        self.logger.info('  Read %d (%.2f Mbp) sequences.' % (len(allSeqs), float(totalBases) / 1e6))

        # write all unbinned sequences
        self.logger.info('Identifying unbinned sequences >= %d bp.' % minSeqLen)
        seqOut = open(outSeqFile, 'w')

        statsOut = open(outStatsFile, 'w')
        statsOut.write('Sequence Id\tLength\tGC\n')

        unbinnedCount = 0
        unbinnedBases = 0
        for seqId, seq in allSeqs.items():
            if seqId not in binnedSeqs:
                if len(seq) >= minSeqLen:
                    unbinnedCount += 1
                    seqOut.write('>' + seqId + '\n')
                    seqOut.write(seq + '\n')

                    unbinnedBases += len(seq)

                    a, c, g, t = baseCount(seq)

                    statsOut.write('%s\t%d\t%.2f\n' % (seqId, len(seq), float(g + c) * 100 / (a + c + g + t)))

        seqOut.close()
        statsOut.close()

        self.logger.info('  Identified %d (%.2f Mbp) unbinned sequences.' % (unbinnedCount, float(unbinnedBases) / 1e6))

        self.logger.info('Percentage of unbinned sequences: %.2f%%' % (unbinnedCount * 100.0 / len(allSeqs)))
        self.logger.info('Percentage of unbinned bases: %.2f%%' % (unbinnedBases * 100.0 / totalBases))
 class Merger():
    def __init__(self):
        self.logger = logging.getLogger('timestamp')

    def run(self, binFiles, outDir, hmmTableFile,
                binIdToModels, binIdToBinMarkerSets,
                minDeltaComp, maxDeltaCont,
                minMergedComp, maxMergedCont):
        checkDirExists(outDir)

        self.logger.info('Comparing marker sets between all pairs of bins.')

        # ensure all bins are using the same marker set
        markerGenesI = binIdToBinMarkerSets[list(binIdToBinMarkerSets.keys())[0]].mostSpecificMarkerSet().getMarkerGenes()
        for binIdJ in binIdToBinMarkerSets:
            if markerGenesI != binIdToBinMarkerSets[binIdJ].mostSpecificMarkerSet().getMarkerGenes():
                self.logger.error('All bins must use the same marker set to assess potential mergers.')
                sys.exit(1)

        # parse HMM information
        resultsParser = ResultsParser(binIdToModels)

        # get HMM hits to each bin
        resultsParser.parseBinHits(outDir, hmmTableFile)

        # determine union and intersection of marker sets for each pair of bins
        outputFile = os.path.join(outDir, "merger.tsv")
        fout = open(outputFile, 'w')
        fout.write('Bin Id 1\tBin Id 2')
        fout.write('\tBin 1 completeness\tBin 1 contamination')
        fout.write('\tBin 2 completeness\tBin 2 contamination')
        fout.write('\tDelta completeness\tDelta contamination\tMerger delta')
        fout.write('\tMerged completeness\tMerged contamination\n')

        binMarkerHits = resultsParser.results
        binIds = sorted(binMarkerHits.keys())
        for i in range(0, len(binMarkerHits)):
            binIdI = binIds[i]

            geneCountsI = binMarkerHits[binIdI].geneCounts(binIdToBinMarkerSets[binIdI].mostSpecificMarkerSet(), binMarkerHits[binIdI].markerHits, True)
            completenessI, contaminationI = geneCountsI[6:8]

            for j in range(i + 1, len(binMarkerHits)):
                binIdJ = binIds[j]

                geneCountsJ = binMarkerHits[binIdJ].geneCounts(binIdToBinMarkerSets[binIdJ].mostSpecificMarkerSet(), binMarkerHits[binIdJ].markerHits, True)
                completenessJ, contaminationJ = geneCountsJ[6:8]

                # merge together hits from both bins and calculate completeness and contamination
                mergedHits = {}
                for markerId, hits in binMarkerHits[binIdI].markerHits.items():
                    mergedHits[markerId] = list(hits)

                for markerId, hits in binMarkerHits[binIdJ].markerHits.items():
                    if markerId in mergedHits:
                        mergedHits[markerId].extend(hits)
                    else:
                        mergedHits[markerId] = hits

                geneCountsMerged = binMarkerHits[binIdI].geneCounts(binIdToBinMarkerSets[binIdJ].mostSpecificMarkerSet(), mergedHits, True)
                completenessMerged, contaminationMerged = geneCountsMerged[6:8]

                if not (completenessMerged >= minMergedComp and contaminationMerged < maxMergedCont):
                    continue

                # calculate merged statistics
                deltaComp = completenessMerged - max(completenessI, completenessJ)
                deltaCont = contaminationMerged - max(contaminationI, contaminationJ)
                delta = deltaComp - deltaCont

                if deltaComp >= minDeltaComp and deltaCont < maxDeltaCont:
                    fout.write('%s\t%s\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n' %
                                                                        (binIdI, binIdJ,
                                                                         completenessI, contaminationI,
                                                                         completenessJ, contaminationJ,
                                                                         deltaComp, deltaCont, delta,
                                                                         completenessMerged, contaminationMerged))

        fout.close()

        return outputFile
 class Profile():
    def __init__(self):
        self.logger = logging.getLogger('timestamp')

    def run(self, coverageFile, outFile, bTabTable):
        checkFileExists(coverageFile)

        # get number of reads mapped to each bin
        self.logger.info('Determining number of reads mapped to each bin.')

        readsMappedToBin = {}
        binSize = {}
        totalMappedReads = {}
        bHeader = True
        for line in open(coverageFile):
            if bHeader:
                bHeader = False
                continue

            lineSplit = line.split('\t')

            # seqId = lineSplit[0]
            binId = lineSplit[1]

            seqLen = int(lineSplit[2])
            binSize[binId] = binSize.get(binId, 0) + seqLen

            if binId not in readsMappedToBin:
                readsMappedToBin[binId] = {}

            for i in range(3, len(lineSplit), 3):
                bamId = lineSplit[i]
                mappedReads = int(lineSplit[i + 2])

                totalMappedReads[bamId] = totalMappedReads.get(bamId, 0) + mappedReads
                readsMappedToBin[binId][bamId] = readsMappedToBin[binId].get(bamId, 0) + mappedReads

        # calculate percentage of mapped reads to binned populations
        perMappedReads = {}
        normBinCoverage = {}
        sumNormBinCoverage = {}
        for binId, bamIds in readsMappedToBin.items():
            perMappedReads[binId] = {}
            normBinCoverage[binId] = {}

            for bamId in bamIds:
                perMR = float(readsMappedToBin[binId][bamId]) / totalMappedReads[bamId]
                perMappedReads[binId][bamId] = perMR

                if binId == DefaultValues.UNBINNED:
                    continue

                normCoverage = perMR / binSize[binId]
                normBinCoverage[binId][bamId] = normCoverage
                sumNormBinCoverage[bamId] = sumNormBinCoverage.get(bamId, 0) + normCoverage

        for binId, bamIds in normBinCoverage.items():
            for bamId in bamIds:
                if sumNormBinCoverage[bamId] != 0:
                    normBinCoverage[binId][bamId] /= sumNormBinCoverage[bamId]
                else:
                    normBinCoverage[binId][bamId] = 0

        # write community profile
        oldStdOut = reassignStdOut(outFile)

        sortedBinIds = sorted(readsMappedToBin.keys())
        sortedBamIds = sorted(readsMappedToBin[sortedBinIds[0]].keys())

        header = ['Bin Id', 'Bin size (Mbp)']
        for bamId in sortedBamIds:
            header += [bamId + ': mapped reads']
            header += [bamId + ': % mapped reads']
            header += [bamId + ': % binned populations']
            header += [bamId + ': % community']

        if bTabTable:
            print('\t'.join(header))
        else:
            pTable = prettytable.PrettyTable(header)
            pTable.float_format = '.2'
            pTable.align = 'c'
            pTable.align[header[0]] = 'l'
            pTable.hrules = prettytable.FRAME
            pTable.vrules = prettytable.NONE

        for binId in sortedBinIds:
            row = [binId]
            row += [float(binSize[binId]) / 1e6]

            for bamId in sortedBamIds:
                row += [readsMappedToBin[binId][bamId]]
                row += [perMappedReads[binId][bamId] * 100.0]

                if DefaultValues.UNBINNED in perMappedReads:
                    unbinnedPercentage = perMappedReads[DefaultValues.UNBINNED][bamId]
                else:
                    unbinnedPercentage = 0

                if binId == DefaultValues.UNBINNED:
                    row += ['NA']
                    row += [unbinnedPercentage * 100.0]
                else:
                    row += [normBinCoverage[binId][bamId] * 100.0]
                    row += [normBinCoverage[binId][bamId] * 100.0 * (1.0 - unbinnedPercentage)]

            if bTabTable:
                print('\t'.join(list(map(str, row))))
            else:
                pTable.add_row(row)

        if not bTabTable:
            print(pTable.get_string())

        restoreStdOut(outFile, oldStdOut)
 class BinTools():
    """Functions for exploring and modifying bins."""
    def __init__(self, threads=1):
        self.logger = logging.getLogger('timestamp')

    def __removeSeqs(self, seqs, seqsToRemove):
        """Remove sequences. """
        missingSeqIds = set(seqsToRemove).difference(set(seqs.keys()))
        if len(missingSeqIds) > 0:
            self.logger.error('Missing sequence(s) specified for removal: ' + ', '.join(missingSeqIds) + '\n')
            sys.exit(1)

        for seqId in seqsToRemove:
            seqs.pop(seqId)

    def __addSeqs(self, seqs, refSeqs, seqsToAdd):
        """Add sequences. """
        missingSeqIds = set(seqsToAdd).difference(set(refSeqs.keys()))
        if len(missingSeqIds) > 0:
            self.logger.error('Missing sequence(s) specified for addition: ' + ', '.join(missingSeqIds) + '\n')
            sys.exit(1)

        for seqId in seqsToAdd:
            seqs[seqId] = refSeqs[seqId]

    def modify(self, binFile, seqFile, seqsToAdd, seqsToRemove, outputFile):
        """Add and remove sequences from a file."""
        binSeqs = readFasta(binFile)

        # add sequences to bin
        if seqsToAdd != None:
            refSeqs = readFasta(seqFile)
            self.__addSeqs(binSeqs, refSeqs, seqsToAdd)

        # remove sequences from bin
        if seqsToRemove != None:
            self.__removeSeqs(binSeqs, seqsToRemove)

        # save modified bin
        writeFasta(binSeqs, outputFile)

    def removeOutliers(self, binFile, outlierFile, outputFile):
        """Remove sequences specified as outliers in the provided file."""

        binSeqs = readFasta(binFile)
        binIdToModify = binIdFromFilename(binFile)

        # get files to remove
        checkFileExists(outlierFile)
        seqsToRemove = []
        bHeader = True
        for line in open(outlierFile):
            if bHeader:
                bHeader = False
                continue

            lineSplit = line.split('\t')
            binId = lineSplit[0]

            if binId == binIdToModify:
                seqId = lineSplit[1]
                seqsToRemove.append(seqId)

        # remove sequences from bin
        if len(seqsToRemove) > 0:
            self.__removeSeqs(binSeqs, seqsToRemove)

        # save modified bin
        writeFasta(binSeqs, outputFile)

    def unique(self, binFiles):
        """Check if sequences are assigned to multiple bins."""

        # read sequence IDs from all bins,
        # while checking for duplicate sequences within a bin
        binSeqs = {}
        for f in binFiles:
            binId = binIdFromFilename(f)

            if f.endswith('.gz'):
                openFile = gzip.open
            else:
                openFile = open

            seqIds = set()
            for line in openFile(f):
                if line[0] == '>':
                    seqId = line[1:].split(None, 1)[0]

                    if seqId in seqIds:
                        print ('  [Warning] Sequence %s found multiple times in bin %s.' % (seqId, binId))
                    seqIds.add(seqId)

            binSeqs[binId] = seqIds

        # check for sequences assigned to multiple bins
        bDuplicates = False
        binIds = list(binSeqs.keys())
        for i in range(0, len(binIds)):
            for j in range(i + 1, len(binIds)):
                seqInter = set(binSeqs[binIds[i]]).intersection(set(binSeqs[binIds[j]]))

                if len(seqInter) > 0:
                    bDuplicates = True
                    print ('  Sequences shared between %s and %s: ' % (binIds[i], binIds[j]))
                    for seqId in seqInter:
                        print ('    ' + seqId)
                    print ('')

        if not bDuplicates:
            print ('  No sequences assigned to multiple bins.')

    def gcDist(self, seqs):
        """GC statistics for bin."""
        GCs = []
        gcTotal = 0
        basesTotal = 0
        for _, seq in seqs.items():
            a, c, g, t = baseCount(seq)
            gc = g + c
            bases = a + c + g + t

            GCs.append(float(gc) / (bases))

            gcTotal += gc
            basesTotal += bases

        meanGC = float(gcTotal) / basesTotal
        deltaGCs = np.array(GCs) - meanGC

        return meanGC, deltaGCs, GCs

    def codingDensityDist(self, seqs, prodigalParser):
        """Coding density statistics for bin."""
        CDs = []

        codingBasesTotal = 0
        basesTotal = 0
        for seqId, seq in seqs.items():
            codingBases = prodigalParser.codingBases(seqId)

            CDs.append(float(codingBases) / len(seq))
            codingBasesTotal += codingBases
            basesTotal += len(seq)

        meanCD = float(codingBasesTotal) / basesTotal
        deltaCDs = np.array(CDs) - meanCD

        return meanCD, deltaCDs, CDs

    def binTetraSig(self, seqs, tetraSigs):
        """Tetranucleotide signature for bin. """
        binSize = 0
        for _, seq in seqs.items():
            binSize += len(seq)

        bInit = True
        for seqId, seq in seqs.items():
            weightedTetraSig = tetraSigs[seqId] * (float(len(seq)) / binSize)
            if bInit:
                binSig = weightedTetraSig
                bInit = False
            else:
                binSig += weightedTetraSig

        return binSig

    def tetraDiffDist(self, seqs, genomicSig, tetraSigs, binSig):
        """TD statistics for bin."""
        deltaTDs = np.zeros(len(seqs))
        for i, seqId in enumerate(seqs.keys()):
            deltaTDs[i] = genomicSig.distance(tetraSigs[seqId], binSig)

        return np.mean(deltaTDs), deltaTDs

    def identifyOutliers(self, outDir, binFiles, tetraProfileFile, distribution, reportType, outputFile):
        """Identify sequences that are outliers."""

        self.logger.info('Reading reference distributions.')
        gcBounds = readDistribution('gc_dist')
        cdBounds = readDistribution('cd_dist')
        tdBounds = readDistribution('td_dist')

        fout = open(outputFile, 'w')
        fout.write('Bin Id\tSequence Id\tSequence length\tOutlying distributions')
        fout.write('\tSequence GC\tMean bin GC\tLower GC bound (%s%%)\tUpper GC bound (%s%%)' % (distribution, distribution))
        fout.write('\tSequence CD\tMean bin CD\tLower CD bound (%s%%)' % distribution)
        fout.write('\tSequence TD\tMean bin TD\tUpper TD bound (%s%%)\n' % distribution)

        processedBins = 0
        for binFile in binFiles:
            binId = binIdFromFilename(binFile)

            processedBins += 1
            self.logger.info('Finding outliers in %s (%d of %d).' % (binId, processedBins, len(binFiles)))

            seqs = readFasta(binFile)

            meanGC, deltaGCs, seqGC = self.gcDist(seqs)

            genomicSig = GenomicSignatures(K=4, threads=1)
            tetraSigs = genomicSig.read(tetraProfileFile)
            binSig = self.binTetraSig(seqs, tetraSigs)
            meanTD, deltaTDs = self.tetraDiffDist(seqs, genomicSig, tetraSigs, binSig)

            gffFile = os.path.join(outDir, 'bins', binId, DefaultValues.PRODIGAL_GFF)
            if not os.path.exists(gffFile):
                self.logger.error('Missing gene feature file (%s). This plot if not compatible with the --genes option.\n' % DefaultValues.PRODIGAL_GFF)
                sys.exit(1)

            prodigalParser = ProdigalGeneFeatureParser(gffFile)
            meanCD, deltaCDs, CDs = self.codingDensityDist(seqs, prodigalParser)

            # find keys into GC and CD distributions
            closestGC = findNearest(np.array(list(gcBounds.keys())), meanGC)
            sampleSeqLen = list(gcBounds[closestGC].keys())[0]
            d = gcBounds[closestGC][sampleSeqLen]
            gcLowerBoundKey = findNearest(list(d.keys()), (100 - distribution) / 2.0)
            gcUpperBoundKey = findNearest(list(d.keys()), (100 + distribution) / 2.0)

            closestCD = findNearest(np.array(list(cdBounds.keys())), meanCD)
            sampleSeqLen = list(cdBounds[closestCD].keys())[0]
            d = cdBounds[closestCD][sampleSeqLen]
            cdLowerBoundKey = findNearest(list(d.keys()), (100 - distribution) / 2.0)

            tdBoundKey = findNearest(list(tdBounds[list(tdBounds.keys())[0]].keys()), distribution)

            index = 0
            for seqId, seq in seqs.items():
                seqLen = len(seq)

                # find GC, CD, and TD bounds
                closestSeqLen = findNearest(list(gcBounds[closestGC].keys()), seqLen)
                gcLowerBound = gcBounds[closestGC][closestSeqLen][gcLowerBoundKey]
                gcUpperBound = gcBounds[closestGC][closestSeqLen][gcUpperBoundKey]

                closestSeqLen = findNearest(list(cdBounds[closestCD].keys()), seqLen)
                cdLowerBound = cdBounds[closestCD][closestSeqLen][cdLowerBoundKey]

                closestSeqLen = findNearest(list(tdBounds.keys()), seqLen)
                tdBound = tdBounds[closestSeqLen][tdBoundKey]

                outlyingDists = []
                if deltaGCs[index] < gcLowerBound or deltaGCs[index] > gcUpperBound:
                    outlyingDists.append('GC')

                if deltaCDs[index] < cdLowerBound:
                    outlyingDists.append('CD')

                if deltaTDs[index] > tdBound:
                    outlyingDists.append('TD')

                if (reportType == 'any' and len(outlyingDists) >= 1) or (reportType == 'all' and len(outlyingDists) == 3):
                    fout.write(binId + '\t' + seqId + '\t%d' % len(seq) + '\t' + ','.join(outlyingDists))
                    fout.write('\t%.1f\t%.1f\t%.1f\t%.1f' % (seqGC[index] * 100, meanGC * 100, (meanGC + gcLowerBound) * 100, (meanGC + gcUpperBound) * 100))
                    fout.write('\t%.1f\t%.1f\t%.1f' % (CDs[index] * 100, meanCD * 100, (meanCD + cdLowerBound) * 100))
                    fout.write('\t%.3f\t%.3f\t%.3f' % (deltaTDs[index], meanTD, tdBound) + '\n')

                index += 1

        fout.close() class SSU_Finder(object):
    def __init__(self, threads):
        self.logger = logging.getLogger('timestamp')

        self.totalThreads = threads

        self.bacteriaModelFile = os.path.join(DefaultValues.CHECKM_DATA_DIR, 'hmms_ssu', 'SSU_bacteria.hmm')
        self.archaeaModelFile = os.path.join(DefaultValues.CHECKM_DATA_DIR, 'hmms_ssu', 'SSU_archaea.hmm')
        self.eukModelFile = os.path.join(DefaultValues.CHECKM_DATA_DIR, 'hmms_ssu', 'SSU_euk.hmm')

    def __hmmSearch(self, seqFile, evalue, outputPrefix):
        if seqFile.endswith('gz'):
            pipe = 'zcat ' + seqFile + ' | '
        else:
            pipe = 'cat ' + seqFile + ' | '

        self.logger.info('Identifying bacterial 16S.')
        os.system(pipe + 'nhmmer --noali --cpu ' + str(self.totalThreads) + ' -o ' + outputPrefix + '.bacteria.txt --tblout ' + outputPrefix + '.bacteria.table.txt -E ' + str(evalue) + ' ' + self.bacteriaModelFile + ' -')

        self.logger.info('Identifying archaeal 16S.')
        os.system(pipe + 'nhmmer --noali --cpu ' + str(self.totalThreads) + ' -o ' + outputPrefix + '.archaea.txt --tblout ' + outputPrefix + '.archaea.table.txt -E ' + str(evalue) + ' ' + self.archaeaModelFile + ' -')

        self.logger.info('Identifying eukaryotic 18S.')
        os.system(pipe + 'nhmmer --noali --cpu ' + str(self.totalThreads) + ' -o ' + outputPrefix + '.euk.txt --tblout ' + outputPrefix + '.euk.table.txt -E ' + str(evalue) + ' ' + self.eukModelFile + ' -')

    def __readHits(self, resultsFile, domain, evalueThreshold):
        """Parse hits from nhmmer output."""
        seqInfo = {}

        bReadHit = False
        for line in open(resultsFile):
            if line[0:2] == '>>':
                lineSplit = line.split()
                seqId = lineSplit[1]
                bReadHit = True
                hitCounter = 0
                continue
            elif line.strip() == '':
                bReadHit = False
            elif bReadHit:
                hitCounter += 1
                if hitCounter >= 3:
                    lineSplit = line.split()

                    iEvalue = lineSplit[3]
                    aliFrom = int(lineSplit[7])
                    aliTo = int(lineSplit[8])

                    revComp = False
                    if aliFrom > aliTo:
                        revComp = True
                        aliFrom, aliTo = aliTo, aliFrom

                    alignLen = int(aliTo) - int(aliFrom)

                    if float(iEvalue) <= evalueThreshold:
                        seqInfo[seqId] = seqInfo.get(seqId, []) + [[domain, iEvalue, str(aliFrom), str(aliTo), str(alignLen), str(revComp)]]

        return seqInfo

    def __addHit(self, hits, seqId, info, concatenateThreshold):
        """Add hits from individual HMMs and concatenate nearby hits."""

        # check if this is the first hit to this sequence
        if seqId not in hits:
            hits[seqId] = info
            return

        # check if hits to sequence are close enough to warrant concatenating them,
        # otherwise record both hits
        baseSeqId = seqId
        index = 1
        bConcatenate = False
        concateSeqId = seqId
        while True:
            # if hits overlap then retain only the longest
            startNew = int(info[2])
            endNew = int(info[3])
            bRevNew = bool(info[5])

            start = int(hits[seqId][2])
            end = int(hits[seqId][3])
            bRev = bool(info[5])

            # check if hits should be concatenated
            if abs(start - endNew) < concatenateThreshold and bRevNew == bRev:
                # new hit closely preceded old hit and is on same strand
                del hits[seqId]
                info[2] = str(startNew)
                info[3] = str(end)
                info[4] = str(end - startNew)
                hits[concateSeqId] = info
                bConcatenate = True

            elif abs(startNew - end) < concatenateThreshold and bRevNew == bRev:
                # new hit closely follows old hit and is on same strand
                del hits[seqId]
                info[2] = str(start)
                info[3] = str(endNew)
                info[4] = str(endNew - start)
                hits[concateSeqId] = info
                bConcatenate = True

            index += 1
            newSeqId = baseSeqId + '-#' + str(index)
            if bConcatenate:
                if newSeqId in hits:
                    seqId = newSeqId  # see if other sequences concatenate
                else:
                    break
            else:
                # hits are not close enough to concatenate
                if newSeqId in hits:
                    seqId = newSeqId  # see if the new hit overlaps with this
                    concateSeqId = newSeqId
                else:
                    hits[newSeqId] = info
                    break

    def __addDomainHit(self, hits, baseSeqId, info):
        """Add hits from different domain models and concatenate nearby hits."""
        
        startNew = int(info[2])
        endNew = int(info[3])
        lengthNew = int(info[4])

        # get all 16S hits for this contig
        curSeqIds = set()
        for curSeqId in hits:
            if baseSeqId in curSeqId:
                curSeqIds.add(curSeqId)

        bNewHit = True
        seqToRemove = set()
        for curSeqId in curSeqIds:
            # if hits overlap then retain only the longest
            start = int(hits[curSeqId][2])
            end = int(hits[curSeqId][3])
            length = int(hits[curSeqId][4])

            if (startNew <= start and endNew >= start) or (start <= startNew and end >= startNew):
                if lengthNew > length:
                    seqToRemove.add(curSeqId)
                else:
                    bNewHit = False
                    break

        if bNewHit:
            max_seq_num = 0
            for curSeqId in curSeqIds:
                if curSeqId in seqToRemove:
                    del hits[curSeqId]
                elif '-#' in curSeqId:
                    seq_num = int(curSeqId.split('-#')[1])
                    if seq_num > max_seq_num:
                        max_seq_num = seq_num

            hits['{}-#{}'.format(baseSeqId, max_seq_num+1)] = info

    def run(self, contigFile, binFiles, outputDir, evalueThreshold, concatenateThreshold):
        # make sure output directory exists
        if not os.path.exists(outputDir):
            os.makedirs(outputDir)

        # get bin id of binned contigs
        self.logger.info('Determining bin assignment of sequences.')
        seqIdToBinId = {}
        for f in binFiles:
            binId = binIdFromFilename(f)
            seqIds = readFastaSeqIds(f)
            for seqId in seqIds:
                seqIdToBinId[seqId] = binId

        # identify 16S reads from contigs/scaffolds
        self.logger.info('Identifying SSU rRNAs on sequences.')
        self.__hmmSearch(contigFile, evalueThreshold, os.path.join(outputDir, 'ssu'))

        # read HMM hits
        hitsPerDomain = {}
        for domain in ['archaea', 'bacteria', 'euk']:
            hits = {}

            seqInfo = self.__readHits(os.path.join(outputDir, 'ssu' + '.' + domain + '.txt'), domain, evalueThreshold)
            if len(seqInfo) > 0:
                for seqId, seqHits in seqInfo.items():
                    for hit in seqHits:
                        self.__addHit(hits, seqId, hit, concatenateThreshold)

            hitsPerDomain[domain] = hits

        # find best domain hit for each sequence
        bestHits = {}
        for _, hits in hitsPerDomain.items():
            for seqId, info in hits.items():
                if '-#' in seqId:
                    seqId = seqId[0:seqId.rfind('-#')]

                self.__addDomainHit(bestHits, seqId, info)
                
        # relabel hits
        newBestHits = {}
        num_hits = defaultdict(int)
        for seqId, seqInfo in bestHits.items():
            if '-#' in seqId:
                seqId = seqId[0:seqId.rfind('-#')]
            
            if seqId not in num_hits:
                newBestHits[seqId] = seqInfo
            else:
                hitNum = num_hits[seqId]
                newBestHits['{}-#{}'.format(seqId, hitNum)] = seqInfo
            
            num_hits[seqId] += 1
            
        bestHits = newBestHits

        # write summary file and putative SSU rRNAs to file
        summaryFile = os.path.join(outputDir, 'ssu_summary.tsv')
        summaryOut = open(summaryFile, 'w')
        summaryOut.write('Bin Id\tSeq. Id\tHMM\ti-Evalue\tStart hit\tEnd hit\t16S/18S gene length\tRev. Complement\tSequence length\n')

        seqFile = os.path.join(outputDir, 'ssu.fna')
        seqOut = open(seqFile, 'w')

        seqs = readFasta(contigFile)

        hitsToBins = {}
        for seqId in bestHits:
            origSeqId = seqId
            if '-#' in seqId:
                seqId = seqId[0:seqId.rfind('-#')]

            if seqId in seqIdToBinId:
                binId = seqIdToBinId[seqId]
            else:
                binId = DefaultValues.UNBINNED

            seqInfo = [origSeqId] + bestHits[origSeqId]
            hitsToBins[binId] = hitsToBins.get(binId, []) + [seqInfo]

        for binId in sorted(hitsToBins.keys()):
            for seqInfo in hitsToBins[binId]:
                seqId = seqInfo[0]
                if '-#' in seqId:
                    seqId = seqId[0:seqId.rfind('-#')]

                seq = seqs[seqId]
                summaryOut.write(binId + '\t' + '\t'.join(seqInfo) + '\t' + str(len(seq)) + '\n')
                seqOut.write('>' + binId + DefaultValues.SEQ_CONCAT_CHAR + seqInfo[0] + '\n')
                seqOut.write(seq[int(seqInfo[3]) - 1:int(seqInfo[4]) - 1] + '\n')

        summaryOut.close()
        seqOut.close()

        self.logger.info('Identified ' + str(len(bestHits)) + ' putative SSU genes.')
        self.logger.info('Summary of identified hits written to: ' + summaryFile)
        self.logger.info('SSU sequences written to: ' + seqFile) def makeSurePathExists(path):
    """Create directory if it does not exist."""
    if not path:
        return

    try:
        os.makedirs(path)
    except OSError as exception:
        if exception.errno != errno.EEXIST:
            logger = logging.getLogger('timestamp')
            logger.error('Specified path does not exist: ' + path + '\n')
            sys.exit(1)
 def checkFileExists(inputFile):
    """Check if file exists."""
    if not os.path.exists(inputFile):
        logger = logging.getLogger('timestamp')
        logger.error('Input file does not exists: ' + inputFile + '\n')
        sys.exit(1)
 def binIdFromFilename(filename):
    """Extract bin id from bin filename."""
    binId = os.path.basename(filename)
    binId = os.path.splitext(binId)[0]

    return binId
 def reassignStdOut(outFile):
    """Redirect standard out to a file."""
    oldStdOut = sys.stdout
    if(outFile != ''):
        try:
            # redirect stdout to a file
            sys.stdout = open(outFile, 'w')
        except:
            logger = logging.getLogger('timestamp')
            logger.error("Error diverting stdout to file: " + outFile)
            sys.exit(1)

    return oldStdOut
 def restoreStdOut(outFile, oldStdOut):
    """Redirect standard out back to system standard out."""
    if(outFile != ''):
        try:
            # redirect stdout to a file
            sys.stdout.close()
            sys.stdout = oldStdOut
        except:
            logger = logging.getLogger('timestamp')
            logger.error("Error restoring stdout ", outFile)
            sys.exit(1)
 def getBinIdsFromOutDir(outDir):
    """Get bin ids."""
    binIds = []

    binDir = os.path.join(outDir, 'bins')

    for f in os.listdir(binDir):
        if os.path.isdir(os.path.join(binDir, f)) and f != 'storage':
            binIds.append(f)

    return binIds
 def checkDirExists(inputDir):
    """Check if directory exists."""
    if not os.path.exists(inputDir):
        logger = logging.getLogger('timestamp')
        logger.error('Input directory does not exists: ' + inputDir + '\n')
        sys.exit(1)
 class GcPlots(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)

    def plot(self, fastaFile, distributionsToPlot):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)

        axesHist = self.fig.add_subplot(121)
        axesDeltaGC = self.fig.add_subplot(122)

        self.plotOnAxes(fastaFile, distributionsToPlot, axesHist, axesDeltaGC)

        self.fig.tight_layout(pad=1, w_pad=1)
        self.draw()

    def plotOnAxes(self, fastaFile, distributionsToPlot, axesHist, axesDeltaGC):
        # Read reference distributions from file
        dist = readDistribution('gc_dist')

        # get GC for windows
        seqs = readFasta(fastaFile)

        data = []
        seqLens = []
        for _, seq in seqs.items():
            start = 0
            end = self.options.gc_window_size

            seqLen = len(seq)
            seqLens.append(seqLen)

            while(end < seqLen):
                a, c, g, t = baseCount(seq[start:end])
                try:
                    data.append(float(g + c) / (a + c + g + t))
                except:
                    # it is possible to reach a long stretch of
                    # N's that causes a division by zero error

                    pass

                start = end
                end += self.options.gc_window_size

        if len(data) == 0:
            axesHist.set_xlabel('[Error] No seqs >= %d, the specified window size' % self.options.gc_window_size)
            return

        # Histogram plot
        bins = [0.0]
        binWidth = self.options.gc_bin_width
        binEnd = binWidth
        while binEnd <= 1.0:
            bins.append(binEnd)
            binEnd += binWidth

        axesHist.hist(data, bins=bins, density=True, color=(0.5, 0.5, 0.5))
        axesHist.set_xlabel('% GC')
        axesHist.set_ylabel('% windows (' + str(self.options.gc_window_size) + ' bp)')

        # Prettify plot
        for a in axesHist.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axesHist.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axesHist.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axesHist.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axesHist.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        # get GC bin statistics
        binTools = BinTools()
        meanGC, deltaGCs, _ = binTools.gcDist(seqs)

        # Delta-GC vs Sequence length plot
        axesDeltaGC.scatter(deltaGCs, seqLens, c=abs(deltaGCs), s=10, lw=0.5, cmap='gray_r')
        axesDeltaGC.set_xlabel(r'$\Delta$ GC (mean GC = %.1f%%)' % (meanGC * 100))
        axesDeltaGC.set_ylabel('Sequence length (kbp)')

        _, yMaxSeqs = axesDeltaGC.get_ylim()
        xMinSeqs, xMaxSeqs = axesDeltaGC.get_xlim()

        # plot reference distributions
        for distToPlot in distributionsToPlot:
            closestGC = findNearest(np.array(list(dist.keys())), meanGC)

            # find closest distribution values
            sampleSeqLen = list(dist[closestGC].keys())[0]
            d = dist[closestGC][sampleSeqLen]
            gcLowerBoundKey = findNearest(list(d.keys()), (100 - distToPlot) / 2.0)
            gcUpperBoundKey = findNearest(list(d.keys()), (100 + distToPlot) / 2.0)

            xL = []
            xU = []
            y = []
            for windowSize in dist[closestGC]:
                xL.append(dist[closestGC][windowSize][gcLowerBoundKey])
                xU.append(dist[closestGC][windowSize][gcUpperBoundKey])
                y.append(windowSize)

            # sort by y-values
            sortIndexY = np.argsort(y)
            xL = np.array(xL)[sortIndexY]
            xU = np.array(xU)[sortIndexY]
            y = np.array(y)[sortIndexY]
            axesDeltaGC.plot(xL, y, 'r--', lw=0.5, zorder=0)
            axesDeltaGC.plot(xU, y, 'r--', lw=0.5, zorder=0)

        # ensure y-axis include zero and covers all sequences
        axesDeltaGC.set_ylim([0, yMaxSeqs])

        # ensure x-axis is set appropriately for sequences
        axesDeltaGC.set_xlim([xMinSeqs, xMaxSeqs])

        # draw vertical line at x=0
        axesDeltaGC.vlines(0, 0, yMaxSeqs, linestyle='dashed', color=self.axesColour, zorder=0)

        # Change sequence lengths from bp to kbp
        yticks = axesDeltaGC.get_yticks()
        kbpLabels = []
        for seqLen in yticks:
            label = '%.1f' % (float(seqLen) / 1000)
            label = label.replace('.0', '')  # remove trailing zero
            kbpLabels.append(label)
        axesDeltaGC.set_yticklabels(kbpLabels)

        # Prettify plot
        for a in axesDeltaGC.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axesDeltaGC.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axesDeltaGC.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axesDeltaGC.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axesDeltaGC.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)
 class CodingDensityPlots(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)
        
        self.logger = logging.getLogger('timestamp')

    def plot(self, fastaFile, distributionsToPlot):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)

        axesHist = self.fig.add_subplot(121)
        axesDeltaCD = self.fig.add_subplot(122)

        self.plotOnAxes(fastaFile, distributionsToPlot, axesHist, axesDeltaCD)

        self.fig.tight_layout(pad=1, w_pad=1)
        self.draw()

    def plotOnAxes(self, fastaFile, distributionsToPlot, axesHist, axesDeltaCD):
        # parse Prodigal output
        gffFile = os.path.join(self.options.results_dir, 'bins', binIdFromFilename(fastaFile), DefaultValues.PRODIGAL_GFF)
        if not os.path.exists(gffFile):
            self.logger.error('Missing gene feature file (%s). This plot if not compatible with the --genes option.' % DefaultValues.PRODIGAL_GFF)
            sys.exit(1)

        prodigalParser = ProdigalGeneFeatureParser(gffFile)

        # Read reference distributions from file
        dist = readDistribution('cd_dist')

        # get coding density for windows
        seqs = readFasta(fastaFile)

        data = []
        seqLens = []
        for seqId, seq in seqs.items():
            start = 0
            end = self.options.cd_window_size

            seqLen = len(seq)
            seqLens.append(seqLen)

            while(end < seqLen):
                codingBases = prodigalParser.codingBases(seqId, start, end)

                a, c, g, t = baseCount(seq[start:end])
                data.append(float(codingBases) / (a + c + g + t))

                start = end
                end += self.options.cd_window_size

        if len(data) == 0:
            axesHist.set_xlabel('[Error] No seqs >= %d, the specified window size' % self.options.cd_window_size)
            return

        # Histogram plot
        bins = [0.0]
        binWidth = self.options.cd_bin_width
        binEnd = binWidth
        while binEnd <= 1.0:
            bins.append(binEnd)
            binEnd += binWidth

        axesHist.hist(data, bins=bins, density=True, color=(0.5, 0.5, 0.5))
        axesHist.set_xlabel('% coding density')
        axesHist.set_ylabel('% windows (' + str(self.options.cd_window_size) + ' bp)')

        # Prettify plot
        for a in axesHist.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axesHist.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axesHist.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axesHist.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axesHist.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        # get CD bin statistics
        binTools = BinTools()
        meanCD, deltaCDs, _ = binTools.codingDensityDist(seqs, prodigalParser)

        # Delta-CD vs sequence length plot
        axesDeltaCD.scatter(deltaCDs, seqLens, c=abs(deltaCDs), s=10, lw=0.5, cmap='gray_r')
        axesDeltaCD.set_xlabel(r'$\Delta$ CD (mean coding density = %.1f%%)' % (meanCD * 100))
        axesDeltaCD.set_ylabel('Sequence length (kbp)')

        _, yMaxSeqs = axesDeltaCD.get_ylim()
        xMinSeqs, xMaxSeqs = axesDeltaCD.get_xlim()

        # plot reference distributions
        for distToPlot in distributionsToPlot:
            closestCD = findNearest(np.array(list(dist.keys())), meanCD)

            # find closest distribution values
            sampleSeqLen = list(dist[closestCD].keys())[0]
            d = dist[closestCD][sampleSeqLen]
            cdLowerBoundKey = findNearest(list(d.keys()), (100 - distToPlot) / 2.0)
            cdUpperBoundKey = findNearest(list(d.keys()), (100 + distToPlot) / 2.0)

            xL = []
            xU = []
            y = []
            for windowSize in dist[closestCD]:
                xL.append(dist[closestCD][windowSize][cdLowerBoundKey])
                xU.append(dist[closestCD][windowSize][cdUpperBoundKey])
                y.append(windowSize)

            # sort by y-values
            sortIndexY = np.argsort(y)
            xL = np.array(xL)[sortIndexY]
            xU = np.array(xU)[sortIndexY]
            y = np.array(y)[sortIndexY]
            axesDeltaCD.plot(xL, y, 'r--', lw=0.5, zorder=0)
            axesDeltaCD.plot(xU, y, 'r--', lw=0.5, zorder=0)

        # ensure y-axis include zero and covers all sequences
        axesDeltaCD.set_ylim([0, yMaxSeqs])

        # ensure x-axis is set appropriately for sequences
        axesDeltaCD.set_xlim([xMinSeqs, xMaxSeqs])

        # draw vertical line at x=0
        axesDeltaCD.vlines(0, 0, yMaxSeqs, linestyle='dashed', color=self.axesColour, zorder=0)

        # Change sequence lengths from bp to kbp
        yticks = axesDeltaCD.get_yticks()
        kbpLabels = []
        for seqLen in yticks:
            label = '%.1f' % (float(seqLen) / 1000)
            label = label.replace('.0', '')  # remove trailing zero
            kbpLabels.append(label)
        axesDeltaCD.set_yticklabels(kbpLabels)

        # Prettify plot
        for a in axesDeltaCD.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axesDeltaCD.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axesDeltaCD.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axesDeltaCD.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axesDeltaCD.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)
 class TetraDistPlots(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)

    def plot(self, fastaFile, tetraSigs, distributionsToPlot):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)

        axesHist = self.fig.add_subplot(121)
        axesDeltaTD = self.fig.add_subplot(122)

        self.plotOnAxes(fastaFile, tetraSigs, distributionsToPlot, axesHist, axesDeltaTD)

        self.fig.tight_layout(pad=1, w_pad=1)
        self.draw()

    def plotOnAxes(self, fastaFile, tetraSigs, distributionsToPlot, axesHist, axesDeltaTD):
        # Read reference distributions from file
        dist = readDistribution('td_dist')

        # get tetranucleotide signature for bin
        seqs = readFasta(fastaFile)

        binTools = BinTools()
        binSig = binTools.binTetraSig(seqs, tetraSigs)

        # get tetranucleotide distances for windows
        genomicSig = GenomicSignatures(K=4, threads=1)

        data = []
        seqLens = []
        deltaTDs = []
        for seqId, seq in seqs.items():
            start = 0
            end = self.options.td_window_size

            seqLen = len(seq)
            seqLens.append(seqLen)
            deltaTDs.append(genomicSig.distance(tetraSigs[seqId], binSig))

            while(end < seqLen):
                windowSig = genomicSig.seqSignature(seq[start:end])
                data.append(genomicSig.distance(windowSig, binSig))

                start = end
                end += self.options.td_window_size

        if len(data) == 0:
            axesHist.set_xlabel('[Error] No seqs >= %d, the specified window size' % self.options.td_window_size)
            return

        deltaTDs = np.array(deltaTDs)

        # Histogram plot
        bins = [0.0]
        binWidth = self.options.td_bin_width
        binEnd = binWidth
        while binEnd <= 1.0:
            bins.append(binEnd)
            binEnd += binWidth

        axesHist.hist(data, bins=bins, density=True, color=(0.5, 0.5, 0.5))
        axesHist.set_xlabel(r'$\Delta$ TD')
        axesHist.set_ylabel('% windows (' + str(self.options.td_window_size) + ' bp)')

        # Prettify plot
        for a in axesHist.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axesHist.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axesHist.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axesHist.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axesHist.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        # get CD bin statistics
        meanTD, deltaTDs = binTools.tetraDiffDist(seqs, genomicSig, tetraSigs, binSig)

        # Delta-TD vs Sequence length plot
        axesDeltaTD.scatter(deltaTDs, seqLens, c=abs(deltaTDs), s=10, lw=0.5, cmap='gray_r')
        axesDeltaTD.set_xlabel(r'$\Delta$ TD (mean TD = %.2f)' % meanTD)
        axesDeltaTD.set_ylabel('Sequence length (kbp)')

        _, yMaxSeqs = axesDeltaTD.get_ylim()
        xMinSeqs, xMaxSeqs = axesDeltaTD.get_xlim()

        # plot reference distributions
        for distToPlot in distributionsToPlot:
            boundKey = findNearest(list(dist[list(dist.keys())[0]].keys()), distToPlot)

            x = []
            y = []
            for windowSize in dist:
                x.append(dist[windowSize][boundKey])
                y.append(windowSize)

            # sort by y-values
            sortIndexY = np.argsort(y)
            x = np.array(x)[sortIndexY]
            y = np.array(y)[sortIndexY]

            # make sure x-values are strictly decreasing as y increases
            # as this is conservative and visually satisfying
            for i in range(0, len(x) - 1):
                for j in range(i + 1, len(x)):
                    if x[j] > x[i]:
                        if j == len(x) - 1:
                            x[j] = x[i]
                        else:
                            x[j] = (x[j - 1] + x[j + 1]) / 2  # interpolate values from neighbours

                        if x[j] > x[i]:
                            x[j] = x[i]

            axesDeltaTD.plot(x, y, 'r--', lw=0.5, zorder=0)

        # ensure y-axis include zero and covers all sequences
        axesDeltaTD.set_ylim([0, yMaxSeqs])

        # ensure x-axis is set appropriately for sequences
        axesDeltaTD.set_xlim([xMinSeqs, xMaxSeqs])

        # draw vertical line at x=0
        axesDeltaTD.vlines(0, 0, yMaxSeqs, linestyle='dashed', color=self.axesColour, zorder=0)

        # Change sequence lengths from bp to kbp
        yticks = axesDeltaTD.get_yticks()
        kbpLabels = []
        for seqLen in yticks:
            label = '%.1f' % (float(seqLen) / 1000)
            label = label.replace('.0', '')  # remove trailing zero
            kbpLabels.append(label)
        axesDeltaTD.set_yticklabels(kbpLabels)

        # Prettify plot
        for a in axesDeltaTD.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axesDeltaTD.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axesDeltaTD.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axesDeltaTD.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axesDeltaTD.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)
 class DistributionPlots(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)
        self.options = options

    def plot(self, fastaFile, tetraSigs, distributionsToPlot):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)

        axesHistGC = self.fig.add_subplot(321)
        axesDeltaGC = self.fig.add_subplot(322)
        axesHistTD = self.fig.add_subplot(323)
        axesDeltaTD = self.fig.add_subplot(324)
        axesHistCD = self.fig.add_subplot(325)
        axesDeltaCD = self.fig.add_subplot(326)

        gcPlots = GcPlots(self.options)
        gcPlots.plotOnAxes(fastaFile, distributionsToPlot, axesHistGC, axesDeltaGC)

        tetraDistPlots = TetraDistPlots(self.options)
        tetraDistPlots.plotOnAxes(fastaFile, tetraSigs, distributionsToPlot, axesHistTD, axesDeltaTD)

        codingDensityPlots = CodingDensityPlots(self.options)
        codingDensityPlots.plotOnAxes(fastaFile, distributionsToPlot, axesHistCD, axesDeltaCD)

        self.fig.tight_layout(pad=1, w_pad=2, h_pad=2)
        self.draw()
 class NxPlot(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)

    def calculateNx(self, x, seqs):
        seqLens = []
        for _, seq in seqs.items():
            seqLens.append(len(seq))

        sumSeqLens = sum(seqLens)

        seqLens.sort(reverse=True)
        x.sort()

        testSum = 0
        xIndex = 0
        nxThreshold = x[xIndex] * sumSeqLens
        nx = []
        for seqLen in seqLens:
            testSum += seqLen

            while(testSum >= nxThreshold):
                nx.append(seqLen)
                if len(nx) == len(x):
                    break

                xIndex += 1
                nxThreshold = x[xIndex] * sumSeqLens

        return nx

    def plot(self, fastaFile):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)
        axes = self.fig.add_subplot(111)

        # calculate Nx
        seqs = readFasta(fastaFile)
        x = np.arange(0, 1.0 + 0.5 * self.options.step_size, self.options.step_size)
        nx = self.calculateNx(x, seqs)

        # Create plot
        axes.plot(x, nx, 'ko-', ms=4)
        axes.set_xlabel('Nx')
        axes.set_ylabel('Sequence length (kbp)')

        # Change sequence lengths from bp to kbp
        yticks = axes.get_yticks()
        kbpLabels = []
        for seqLen in yticks:
            label = '%.1f' % (float(seqLen) / 1000)
            label = label.replace('.0', '')  # remove trailing zero
            kbpLabels.append(label)
        axes.set_yticklabels(kbpLabels)

        # Prettify plot
        for a in axes.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axes.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axes.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axes.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        self.fig.tight_layout(pad=1)
        self.draw()
 class LengthHistogram(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)

    def plot(self, fastaFile):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)
        axes = self.fig.add_subplot(111)

        # calculate sequence lengths (in kb)
        seqs = readFasta(fastaFile)

        seqLens = []
        for seq in seqs.values():
            seqLens.append(float(len(seq)) / 1e3)

        # set unequal bin sizes (in kb)
        bins = [0, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1e12]
        counts, _edges = np.histogram(seqLens, bins=bins)

        # create histogram
        axes.bar(x=np.arange(0.1, len(counts)), height=counts, width=0.8, color=(0.5, 0.5, 0.5))
        axes.set_xlabel('Sequence length (kbp)')
        axes.set_ylabel('Number sequences (out of %d)' % len(seqs))

        # ensure y-axis include zero
        _, end = axes.get_ylim()
        axes.set_ylim([0, end])
        axes.get_yaxis().set_major_locator(MaxNLocator(integer=True))

        # Change sequence lengths from bp to kbp
        axes.set_xlim([0, len(counts)])
        axes.set_xticks(np.arange(0.5, len(counts)))
        axes.set_xticklabels(['<1', '1-2', '2-5', '5-10', '10-20', '20-50', '50-100', '100-200', '200-500', '>500'])

        # Prettify plot
        for a in axes.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in axes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axes.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axes.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in axes.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        self.fig.tight_layout(pad=1)
        self.draw()
 class MarkerGenePosPlot(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)
        
        self.logger = logging.getLogger('timestamp')

    def getMarkerGenesPerSeq(self, markerGeneStats):
        markerGenePos = {}
        markerGeneNum = {}

        for geneId in markerGeneStats:
            scaffoldId = geneId[0:geneId.rfind('_')]

            for markerId, hitList in markerGeneStats[geneId].items():
                for hit in hitList:
                    start = hit[0]
                    end = hit[1]
                markerGenePos[scaffoldId] = markerGenePos.get(scaffoldId, []) + [[geneId, markerId, start, end]]
                markerGeneNum[markerId] = len(hitList)

        return markerGenePos, markerGeneNum

    def roundUpToNearest100(self, x):
        return int(math.ceil(x / 100.0)) * 100

    def plot(self, binFile, markerGeneStats, binStats):
        binId = binIdFromFilename(binFile)

        markerGenesPerSeq, _markerGeneNum = self.getMarkerGenesPerSeq(markerGeneStats)

        if len(markerGenesPerSeq) == 0:
            return False

        # Get length of sequences with one or more marker genes
        seqs = readFasta(binFile)
        seqLens = {}
        longestSeq = 0
        binSize = 0
        for seqId, seq in seqs.items():
            seqLen = len(seq)
            binSize += seqLen

            if seqId not in markerGenesPerSeq:
                continue

            seqLens[seqId] = seqLen
            if seqLen > longestSeq:
                longestSeq = seqLen

        sortedSeqLens = sorted(seqLens.items(), key=operator.itemgetter(1), reverse=True)
        
        result_str = 'Markers reside on {:,} of {:,} sequences which span {:.2f} of {:.2f} ({:.1f}%) Mb'.format(
                            len(seqLens),
                            len(seqs),
                            sum([s for s in seqLens.values()])/1e6,
                            binSize/1e6,
                            sum([s for s in seqLens.values()])*100.0/binSize)
        self.logger.info(result_str)

        MAX_BINS = 100
        plotBinSize = self.roundUpToNearest100(float(longestSeq) / MAX_BINS)
        yLabels = [x[0] for x in sortedSeqLens]

        # get position of genes in bin
        prodigalFastaParser = ProdigalFastaParser()
        geneFile = os.path.join(self.options.results_dir, 'bins', binId, DefaultValues.PRODIGAL_AA)
        genePos = prodigalFastaParser.genePositions(geneFile)

        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)
        yLabelBounds = self.yLabelExtents(yLabels, self.options.font_size)

        heightBottomLabels = 0.4 + self.options.fig_padding  # inches
        widthSideLabel = yLabelBounds.width * self.options.width + self.options.fig_padding  # inches

        widthPerBin = (self.options.width - widthSideLabel - self.options.fig_padding) / MAX_BINS

        titleHeight = 0.2
        HEIGHT_PER_ROW = 0.2
        height = HEIGHT_PER_ROW * len(sortedSeqLens) + heightBottomLabels + self.options.fig_padding + titleHeight
        rowBinHeight = widthPerBin / HEIGHT_PER_ROW

        self.fig.set_size_inches(self.options.width, height)
        axes = self.fig.add_axes([widthSideLabel / self.options.width, heightBottomLabels / height, \
                                                                        1.0 - (widthSideLabel + self.options.fig_padding) / self.options.width, \
                                                                        1.0 - (heightBottomLabels + self.options.fig_padding + titleHeight) / height])

        # set plot axis
        axes.set_xlim([0, MAX_BINS + 0.1])
        axes.set_xlabel('Position ({:,} bp/bin)'.format(plotBinSize))

        axes.set_ylim([0, len(sortedSeqLens)])
        axes.set_yticks(np.arange(0.5, len(sortedSeqLens) + 0.5, 1.0))

        axes.set_yticklabels(yLabels)

        # legend
        colours = [(1.0, 1.0, 1.0), 
                    (127 / 255.0, 201 / 255.0, 127 / 255.0), 
                    (255 / 255.0, 192 / 255.0, 134 / 255.0), 
                    (190 / 255.0, 174 / 255.0, 212 / 255.0), 
                    (0.0, 0.0, 0.0)]
        discreteColourMap = mpl.colors.ListedColormap(colours)
        axisColourMap = self.fig.add_axes([self.options.fig_padding / self.options.width, self.options.fig_padding / height, 0.15, 0.03 * (self.options.width / height)])
        colourBar = mpl.colorbar.ColorbarBase(axisColourMap, cmap=discreteColourMap, norm=mpl.colors.Normalize(vmin=0, vmax=1), orientation='horizontal', drawedges=True)
        colourBar.set_ticks([0.1, 0.3, 0.5, 0.7, 0.9])
        colourBar.set_ticklabels(['0', '1', '2', '3', '4+'])
        colourBar.outline.set_linewidth(0.5)
        colourBar.dividers.set_linewidth(0.5)

        for a in axisColourMap.xaxis.majorTicks:
            a.tick1On = False
            a.tick2On = False

        # plot each bin
        binPosX = 0.5
        for seqId, seqLen in sortedSeqLens:
            markerCount = [0] * int(math.ceil(float(seqLen) / plotBinSize))
            for geneId, _markerGeneId, geneStartPos, _geneEndPos in markerGenesPerSeq[seqId]:
                binPos = int(float(genePos[geneId][0] + geneStartPos) / plotBinSize)
                markerCount[binPos] += 1

            for i in range(0, len(markerCount)):
                if markerCount[i] < len(colours):
                    axes.add_patch(Rectangle((i + 0.1, binPosX - 0.4 * rowBinHeight), 0.8, 0.8 * rowBinHeight, facecolor=colours[markerCount[i]], lw=0.2))
                else:
                    axes.add_patch(Rectangle((i + 0.1, binPosX - 0.4 * rowBinHeight), 0.8, 0.8 * rowBinHeight, facecolor=colours[-1], lw=0.2))

            binPosX += 1.0

        # set plot title
        titleStr = binId
        titleStr += '\n'
        titleStr += '({:.1f}% complete, {:.1f}% contamination)'.format(
                        binStats['Completeness'], 
                        binStats['Contamination'])
        titleStr += '\n'
        titleStr += '({})'.format(result_str)
        axes.set_title(titleStr)

        # Prettify plot
        for a in axes.yaxis.majorTicks:
            a.tick1On = False
            a.tick2On = False

        for a in axes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in axes.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in axes.xaxis.get_ticklines():
            line.set_color(self.axesColour)
            line.set_ms(2)

        for loc, spine in axes.spines.items():
            if loc in ['left', 'right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        self.draw()

        return True
 class GcBiasPlot(AbstractPlot):
    def __init__(self, options):
        AbstractPlot.__init__(self, options)

    def plot(self, binFile, coverageProfile):
        # Set size of figure
        self.fig.clear()
        self.fig.set_size_inches(self.options.width, self.options.height)

        windowAxes = self.fig.add_subplot(121)
        seqAxes = self.fig.add_subplot(122)

        self.plotOnAxes(binFile, coverageProfile, windowAxes, seqAxes)

        self.fig.tight_layout(pad=1)
        self.draw()

    def plotOnAxes(self, binFile, coverageProfile, windowAxes, seqAxes):

        # get GC for windows
        seqs = readFasta(binFile)

        gcProfile = {}
        for seqId, seq in seqs.items():
            start = 0
            end = self.options.window_size

            windowGCs = []
            while(end < len(seq)):
                a, c, g, t = baseCount(seq[start:end])
                windowGCs.append(float(g + c) / (a + c + g + t))

                start = end
                end += self.options.window_size

            a, c, g, t = baseCount(seq)
            seqGC = float(g + c) / (a + c + g + t)
            gcProfile[seqId] = [seqGC, windowGCs]

        # plot GC vs coverage for windows
        gc = []
        coverage = []
        for seqId, gcInfo in gcProfile.items():
            gc += gcInfo[1]
            coverage += coverageProfile[seqId][1]

        windowAxes.scatter(gc, coverage, c=abs(array(coverage)), s=10, lw=0.5, cmap='gray_r')
        windowAxes.set_xlabel('GC (mean = %.1f%%)' % (mean(gc) * 100))
        windowAxes.set_ylabel('Coverage (mean = %.1f)' % mean(coverage))

        # plot linear regression line
        if len(gc) > 1:
            slope, inter = polyfit(gc, coverage, 1)
            fit_fn = poly1d([slope, inter])  # fit_fn is now a function which takes in x and returns an estimate for y
            windowAxes.plot([min(gc), max(gc)], fit_fn([min(gc), max(gc)]), '--r', lw=0.5)
            windowAxes.set_title('GC vs. Coverage\n(window size = %d bp, slope = %.2f)' % (self.options.window_size, slope))
        else:
            # not possible to calculate best fit line
            windowAxes.set_title('GC vs. Coverage\n(window size = %d bp, no best fit line)' % self.options.window_size)

        # Prettify plot
        for a in windowAxes.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in windowAxes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in windowAxes.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in windowAxes.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in windowAxes.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)

        # plot GC vs coverage for entire sequences
        gc = []
        coverage = []
        seqLen = []
        for seqId, gcInfo in gcProfile.items():
            gc.append(gcInfo[0])
            coverage.append(coverageProfile[seqId][0])
            seqLen.append(len(seqs[seqId]))

        # set marker size proportional to sequence length
        markerSize = log(array(seqLen))  # log-scale
        markerSize = (markerSize - min(markerSize)) / max(markerSize)  # normalize between 0 and 1
        markerSize = markerSize * 200 + 10  # normalize between 10 and 200

        seqAxes.scatter(gc, coverage, c=abs(array(coverage)), s=markerSize, lw=0.5, cmap='gray_r')
        seqAxes.set_xlabel('GC (mean = %.1f%%)' % (mean(gc) * 100))
        seqAxes.set_ylabel('Coverage (mean = %.1f)' % mean(coverage))
        seqAxes.set_title('GC vs. Coverage\nIndividual Sequences')

        # Prettify plot
        for a in seqAxes.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in seqAxes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in seqAxes.yaxis.get_ticklines():
            line.set_color(self.axesColour)

        for line in seqAxes.xaxis.get_ticklines():
            line.set_color(self.axesColour)

        for loc, spine in seqAxes.spines.items():
            if loc in ['right', 'top']:
                spine.set_color('none')
            else:
                spine.set_color(self.axesColour)
 def checkNuclotideSeqs(seq_files):
    """Check if files contain sequences in nucleotide space.

    Parameters
    ----------
    seq_files : iterable
        Sequence files to check.

    Returns
    -------
    boolean
        True if files can be treated as containing nucleotide sequences.
    """

    for seq_file in seq_files:
        if os.stat(seq_file).st_size == 0:
            continue
            
        if not isNucleotide(seq_file):
            logger = logging.getLogger('timestamp')
            logger.warning('Expected all files to contain sequences in nucleotide space.')
            logger.warning('File %s appears to contain amino acids sequences.' % seq_file)

    return True
 def checkProteinSeqs(seq_files):
    """Check if files contain sequences in amino acid space.

    Parameters
    ----------
    seq_files : iterable
        Sequence files to check.

    Returns
    -------
    boolean
        True if files can be treated as containing amino acid sequences.
    """

    for seq_file in seq_files:
        if os.stat(seq_file).st_size == 0:
            continue
        
        if isNucleotide(seq_file):
            logger = logging.getLogger('timestamp')
            logger.warning('Expected all files to contain sequences in amino acid space.')
            logger.warning('File %s appears to contain nucleotide sequences.' % seq_file)

    return True
 class DBManager(mm.ManifestManager):

    """Manage all aspects of data location and version control."""
    def __init__(self, set_path=None):
        mm.ManifestManager.__init__(self, timeout=15)
        self.logger = logging.getLogger('timestamp')
        self.config = DBConfig()  # load inbuilt configuration
        self.type = self.config.values["manifestType"]
        
        if set_path:
            self.setRoot(set_path)

        if set_path:
            self.setRoot(set_path)

        # check that the data root is legit
        manifestFile = os.path.join(self.config.values["dataRoot"], mm.__MANIFEST__)
        if not os.path.exists(self.config.values["dataRoot"]) or not os.path.exists(manifestFile):
            self.config.values["dataRoot"] = ""

        if self.config.values["dataRoot"] == "":
            # no data folder set.
            print ("It seems that the CheckM data folder has not been set yet or has been removed. Please run 'checkm data setRoot'.")
            if not self.setRoot():
                print("Sorry, CheckM cannot run without a valid data folder.")

    def runAction(self, action):
        """Main entry point for the updating code"""
        
        if action[0] == "setRoot":
            if len(action) > 1:
                path = self.setRoot(path=action[1])
            else:
                path = self.setRoot()

            if path is None:
                self.logger.info("Data location not changed")
            else:
                self.logger.info("Data location successfully changed to: %s" % path)
        else:
            self.logger.error("Unknown action: %s" % action[0])

    def setRoot(self, path=None):
        """Set the data folder"""
        # check to see we have permission to update the data root
        if not self.config.checkPermissions():
            return None

        # validate the supplied path or prompt for a new one
        path = self.confirmPath(path=path)
        if path is None:
            # The user is not interested
            return None

        # path should be set, exist and be writable
        self.config.values["dataRoot"] = path

        # save the new path
        self.config.setConfig()

        return path

    def confirmPath(self, path=None):
        """Ask the user to supply a path"""
        path_set = False
        minimal = False
        while not path_set:
            if not path:
                path = os.path.join(os.path.expanduser("~"), ".checkm")

            if path.upper() == "ABORT":
                # user has given up
                return None
            else:
                path = os.path.abspath(os.path.expanduser(path))

            print ("")
            if os.path.exists(path):
                # path exists
                if os.access(path, os.W_OK):
                    # path is writable
                    path_set = True
                    print ("Path [%s] exists and you have permission to write to this folder." % path)
                else:
                    print ("Path [%s] exists but you do not have permission to write to this folder." % path)
                    sys.exit(-1)
            else:
                # path does not exist, try to make it
                print ("Path [%s] does not exist so I will attempt to create it" % path)
                try:
                    self.makeSurePathExists(path)
                    print ("Path [%s] has been created and you have permission to write to this folder." % path)
                    path_set = True
                except Exception:
                    print ("Unable to make the folder; error was: %s" % sys.exc_info()[0])
                minimal = True

        # (re)make the manifest file
        print ("(re) creating manifest file (please be patient).")
        self.createManifest(path, self.config.values["localManifestName"])

        return path

    def checkPermissions(self):
        """See if the user has permission to write to the data directory"""
        if not os.access(self.config.values["dataRoot"], os.W_OK):
            print ("You do not seem to have permission to edit the CheckM data folder")
            print ("located at %s" % self.config.values["dataRoot"])
            return False

        return True
 class VerifyEcoli():
    def __init__(self):
        self.logger = logging.getLogger('timestamp')

    def run(self, parser, outputDir):
        """Run standard E. coli genome to verify operation of CheckM."""

        ecoliFile = os.path.join(DefaultValues.CHECKM_DATA_DIR, 'test_data', '637000110.fna')
        checkFileExists(ecoliFile)

        options = Options()
        options.threads = 1
        options.pplacer_threads = 1
        options.extension = 'fna'
        options.bQuiet = True
        options.output_dir = os.path.join(outputDir, 'results')
        if os.path.exists(options.output_dir):
            shutil.rmtree(options.output_dir)
        makeSurePathExists(options.output_dir)

        self.logger.info('[Step 1]: Verifying tree command.')
        options.bKeepAlignment = False
        options.bNucORFs = False
        options.bCalledGenes = False
        options.bReducedTree = True
        options.bin_dir = os.path.join(DefaultValues.CHECKM_DATA_DIR, 'test_data')
        parser.tree(options)
        self.verifyTree(options.output_dir)
        self.logger.info('[Passed]')

        self.logger.info('[Step 2]: Verifying tree_qa command.')
        options.tree_dir = options.output_dir
        options.out_format = 1
        options.file = os.path.join(options.output_dir, 'tree_qa_test.tsv')
        options.bTabTable = True
        parser.treeQA(options)
        self.verifyTreeQA(options.file)
        self.logger.info('[Passed]')

        self.logger.info('[Step 3]: Verifying lineage_set command.')
        options.marker_file = os.path.join(options.output_dir, 'lineage_set_test.tsv')
        options.bForceDomain = False
        options.bootstrap = 0
        options.num_genomes_markers = 30
        options.num_genomes_refine = 5
        options.bNoLineageSpecificRefinement = False

        options.bRequireTaxonomy = False
        options.unique = 10
        options.multi = 10
        parser.lineageSet(options)
        self.verifyLineageSet(options.marker_file, options.bRequireTaxonomy)

        options.bRequireTaxonomy = True
        parser.lineageSet(options)
        self.verifyLineageSet(options.marker_file, options.bRequireTaxonomy)
        self.logger.info('[Passed]')

        self.logger.info('[Step 4]: Verifying analyze command.')
        options.bAlignTopHit = False
        parser.analyze(options)
        self.verifyAnalyze(options.output_dir)
        self.logger.info('[Passed]')

        self.logger.info('[Step 5]: Verifying qa command.')
        options.alignment_file = None
        options.analyze_dir = options.output_dir
        options.out_format = 1
        options.exclude_markers = None
        options.bSkipPseudoGeneCorrection = False
        options.bSkipAdjCorrection = False
        options.file = os.path.join(options.output_dir, 'qa_test.tsv')
        options.bIndividualMarkers = False
        options.bIgnoreThresholds = False
        options.aai_strain = 0.9
        options.e_value = 1e-10
        options.length = 0.7
        options.coverage_file = None
        options.bTabTable = True
        parser.qa(options)
        self.verifyQA(options.file)
        self.logger.info('[Passed]')

    def verifyTree(self, outdir):
        """Verify output of tree command."""

        # verify bin stats using independently verified ground truth values
        binStats = ResultsParser(None).parseBinStats(outdir, DefaultValues.BIN_STATS_PHYLO_OUT)

        np.testing.assert_almost_equal(binStats['637000110']['GC'], 0.508, decimal=3, err_msg="Failed GC test")
        np.testing.assert_almost_equal(binStats['637000110']['GC std'], 0.0, err_msg="Failed GC std test")
        # np.testing.assert_almost_equal(binStats['637000110']['Coding density'], 0.8775, decimal=3, err_msg="Failed coding density test") # depends on exact version of prodigal
        np.testing.assert_almost_equal(binStats['637000110']['# contigs'], 1, err_msg="Failed # contigs test")
        np.testing.assert_almost_equal(binStats['637000110']['# scaffolds'], 1, err_msg="Failed # scaffolds test")
        np.testing.assert_equal(binStats['637000110']['Longest contig'], 4646332, err_msg="Failed longest contig test")
        np.testing.assert_equal(binStats['637000110']['Longest scaffold'], 4646332, err_msg="Failed longest scaffold test")
        # np.testing.assert_equal(binStats['637000110']['# predicted genes'], 4327, err_msg="Failed # predicted genes test") # depends on exact version of prodigal
        np.testing.assert_equal(binStats['637000110']['N50 (contigs)'], 4646332, err_msg="Failed N50 (contigs) test")
        np.testing.assert_equal(binStats['637000110']['N50 (scaffolds)'], 4646332, err_msg="Failed N50 (scaffolds) test")
        np.testing.assert_equal(binStats['637000110']['Genome size'], 4646332, err_msg="Failed genome size test")

        # verify sequence stats using  independently verified ground truth values
        # [The sequence stats file is not generated in CheckM >= v0.9.8 in order to reduce memory requirements.]
        # with open(os.path.join(outdir, 'storage', DefaultValues.SEQ_STATS_PHYLO_OUT), 'r') as f:
        #    s = f.read()
        #    seqStats = ast.literal_eval(s)

        # np.testing.assert_almost_equal(seqStats['637000110']['AC_000091']['GC'], 0.508, decimal=3, err_msg="Failed GC test")
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['Total contig length'], 4646332, err_msg="Failed total contig length test")
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['Coding bases'], 4077069, err_msg="Failed coding bases test") # depends on exact version of prodigal
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['# ORFs'], 4326, err_msg="Failed # genes test") # depends on exact version of prodigal
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['Length'], 4646332, err_msg="Failed length test")
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['# contigs'], 1, err_msg="Failed # contigs test")

    def verifyTreeQA(self, qaTableFile):
        """Verify output of tree QA command."""

        with open(qaTableFile) as f:
            f.readline()  # skip header

            for line in f:
                if line.strip() != '':
                    lineSplit = line.split('\t')

        np.testing.assert_almost_equal(int(lineSplit[0]), 637000110, err_msg="Failed genome ID test")
        np.testing.assert_almost_equal(int(lineSplit[1]), 43, err_msg="Failed # markers test")

        family = None
        if len(lineSplit) >= 4:
            taxonomy = lineSplit[3].split(';')
            if len(taxonomy) >= 5:
                family = lineSplit[3].split(';')[4].strip()
        assert(family == 'f__Enterobacteriaceae')

    def verifyLineageSet(self, markerSetFile, bRequireTaxonomy):
        """Verify output of lineage set command."""

        with open(markerSetFile) as f:
            f.readline()  # skip header

            for line in f:
                if line.strip() != '':
                    lineSplit = line.split('\t')
                    binId = lineSplit[0]
                    _numMarkers = int(lineSplit[1])
                    uid = lineSplit[2]
                    lineage = lineSplit[3]
                    numGenomes = int(lineSplit[4])
                    _markerSet = MarkerSet(uid, lineage, numGenomes, eval(lineSplit[5].rstrip()))

        np.testing.assert_almost_equal(int(binId), 637000110, err_msg="Failed bin ID test")
        if not bRequireTaxonomy:
            # this is unstable as it depends on HMMER and prodigal
            # np.testing.assert_equal(markerSet.numSets(), 336, err_msg="Failed # marker set test")
            # np.testing.assert_equal(markerSet.numMarkers(), 1173, err_msg="Failed # markers test")
            pass
        else:
            # np.testing.assert_equal(markerSet.numSets(), 282, err_msg="Failed # marker set test")
            # np.testing.assert_equal(markerSet.numMarkers(), 1254, err_msg="Failed # markers test")
            pass

    def verifyAnalyze(self, outdir):
        """Verify output of analyze command."""

        # verify bin stats using independently verified ground truth values
        binStats = ResultsParser(None).parseBinStats(outdir, DefaultValues.BIN_STATS_OUT)

        np.testing.assert_almost_equal(binStats['637000110']['GC'], 0.508, decimal=3, err_msg="Failed GC test")
        np.testing.assert_almost_equal(binStats['637000110']['GC std'], 0.0, err_msg="Failed GC std test")
        # np.testing.assert_almost_equal(binStats['637000110']['Coding density'], 0.877, decimal=3, err_msg="Failed coding density test") # depends on exact version of prodigal
        np.testing.assert_equal(binStats['637000110']['# contigs'], 1, err_msg="Failed # contigs test")
        np.testing.assert_equal(binStats['637000110']['# scaffolds'], 1, err_msg="Failed # scaffolds test")
        np.testing.assert_equal(binStats['637000110']['Longest contig'], 4646332, err_msg="Failed longest contig test")
        np.testing.assert_equal(binStats['637000110']['Longest scaffold'], 4646332, err_msg="Failed longest scaffold test")
        # np.testing.assert_equal(binStats['637000110']['# predicted genes'], 4326, err_msg="Failed # predicted genes test") # depends on exact version of prodigal
        np.testing.assert_equal(binStats['637000110']['N50 (contigs)'], 4646332, err_msg="Failed N50 (contigs) test")
        np.testing.assert_equal(binStats['637000110']['N50 (scaffolds)'], 4646332, err_msg="Failed N50 (scaffolds) test")
        np.testing.assert_equal(binStats['637000110']['Genome size'], 4646332, err_msg="Failed genome size test")

        # verify sequence stats using  independently verified ground truth values
        # [The sequence stats file is not generated in CheckM >= v0.9.8 in order to reduce memory requirements.]
        # with open(os.path.join(outdir, 'storage', DefaultValues.SEQ_STATS_OUT), 'r') as f:
        #    s = f.read()
        #    seqStats = ast.literal_eval(s)

        # np.testing.assert_almost_equal(seqStats['637000110']['AC_000091']['GC'], 0.508, decimal=3, err_msg="Failed GC test")
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['Total contig length'], 4646332, err_msg="Failed total contig length test")
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['Coding bases'], 4077069, err_msg="Failed coding bases test") # depends on exact version of prodigal
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['# ORFs'], 4326, err_msg="Failed # genes test") # depends on exact version of prodigal
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['Length'], 4646332, err_msg="Failed length test")
        # np.testing.assert_equal(seqStats['637000110']['AC_000091']['# contigs'], 1, err_msg="Failed # contigs test")

    def verifyQA(self, qaTableFile):
        """Verify output of qa command."""

        with open(qaTableFile) as f:
            f.readline()  # skip header

            for line in f:
                if line.strip() != '':
                    lineSplit = line.split('\t')

        np.testing.assert_equal(int(lineSplit[0]), 637000110, err_msg="Failed genome ID test")
        # np.testing.assert_equal(lineSplit[1], 'f__Enterobacteriaceae', err_msg="Failed lineage test")    # depend on exact version of prodigal
        # np.testing.assert_equal(int(lineSplit[2]), 134, err_msg="Failed # genomes")
        # np.testing.assert_equal(int(lineSplit[3]), 1173, err_msg="Failed # markers test")
        # np.testing.assert_almost_equal(int(lineSplit[4]), 336, err_msg="Failed # marker sets test")
        # np.testing.assert_almost_equal(float(lineSplit[11]), 99.98, decimal=2, err_msg="Failed completeness test")
        # np.testing.assert_almost_equal(float(lineSplit[12]), 0.04, decimal=2, err_msg="Failed contamination test")
import os
import sys
import logging
import numpy as np
from collections import defaultdict
from checkm.defaultValues import DefaultValues
from checkm.timeKeeper import TimeKeeper
from checkm.markerSets import MarkerSetParser, BinMarkerSets
from checkm.resultsParser import ResultsParser
from checkm.hmmerAligner import HmmerAligner
from checkm.markerGeneFinder import MarkerGeneFinder
from checkm.pplacer import PplacerRunner
from checkm.treeParser import TreeParser
from checkm.taxonParser import TaxonParser
from checkm.aminoAcidIdentity import AminoAcidIdentity
from checkm.binStatistics import BinStatistics
from checkm.coverage import Coverage
from checkm.coverageWindows import CoverageWindows
from checkm.genomicSignatures import GenomicSignatures
from checkm.unbinned import Unbinned
from checkm.merger import Merger
from checkm.profile import Profile
from checkm.binTools import BinTools
from checkm.ssuFinder import SSU_Finder
from checkm.common import (makeSurePathExists, 
                            checkFileExists, 
                            binIdFromFilename, 
                            reassignStdOut, 
                            restoreStdOut, 
                            getBinIdsFromOutDir, 
                            checkDirExists)
from checkm.plot.gcPlots import GcPlots
from checkm.plot.codingDensityPlots import CodingDensityPlots
from checkm.plot.tetraDistPlots import TetraDistPlots
from checkm.plot.distributionPlots import DistributionPlots
from checkm.plot.nxPlot import NxPlot
from checkm.plot.lengthHistogram import LengthHistogram
from checkm.plot.markerGenePosPlot import MarkerGenePosPlot
from checkm.plot.gcBiasPlots import GcBiasPlot
from checkm.util.seqUtils import checkNuclotideSeqs, checkProteinSeqs
from checkm.checkmData import DBManager
from checkm.test.test_ecoli import VerifyEcoli

            DBM = DBManager(set_path=action[1])
        else:
            self.logger.error('Path to the CheckM reference data must be specified.')

    def binFiles(self, binFolder, binExtension):
        binFiles = []
        if binFolder is not None:
            all_files = os.listdir(binFolder)
            for f in all_files:
                if f.endswith(binExtension):
                    binFile = os.path.join(binFolder, f)
                    if os.stat(binFile).st_size == 0:
                        self.logger.warning("Skipping bin %s as it has a size of 0 bytes." % f)
                    else:
                        binFiles.append(binFile)

        if not binFiles:
            self.logger.error("No bins found. Check the extension (-x) used to identify bins.")
            sys.exit(1)

        return sorted(binFiles)

    def tree(self, options):
        """Tree command"""
        self.logger.info('[CheckM - tree] Placing bins in reference genome tree.')

        binFiles = self.binFiles(options.bin_dir, options.extension)

        if not options.bCalledGenes:
            if not checkNuclotideSeqs(binFiles):
                return
        else:
            if not checkProteinSeqs(binFiles):
                return

        # setup directory structure
        makeSurePathExists(os.path.join(options.output_dir, 'bins'))
        makeSurePathExists(os.path.join(options.output_dir, 'storage'))

        # find phylogenetically informative genes in genome bins
        mgf = MarkerGeneFinder(options.threads)
        binIdToModels = mgf.find(binFiles,
                                 options.output_dir,
                                 DefaultValues.HMMER_TABLE_PHYLO_OUT,
                                 DefaultValues.HMMER_PHYLO_OUT,
                                 DefaultValues.PHYLO_HMM_MODELS,
                                 options.bKeepAlignment,
                                 options.bNucORFs,
                                 options.bCalledGenes)

        # write model information to file
        markerSetParser = MarkerSetParser(options.threads)
        hmmModelInfoFile = os.path.join(options.output_dir, 'storage', DefaultValues.PHYLO_HMM_MODEL_INFO)
        markerSetParser.writeBinModels(binIdToModels, hmmModelInfoFile)

        # calculate statistics for each genome bin
        
        binStats = BinStatistics(options.threads)
        binStats.calculate(binFiles, options.output_dir, DefaultValues.BIN_STATS_PHYLO_OUT)

        # align identified marker genes
        
        HA = HmmerAligner(options.threads)
        resultsParser = HA.makeAlignmentToPhyloMarkers(options.output_dir,
                                                            DefaultValues.PHYLO_HMM_MODELS,
                                                            DefaultValues.HMMER_TABLE_PHYLO_OUT,
                                                            binIdToModels,
                                                            False,
                                                            DefaultValues.E_VAL,
                                                            DefaultValues.LENGTH,
                                                            False,
                                                            os.path.join(options.output_dir, 'storage', 'tree')
                                                            )

        # place bins into genome tree
        
        pplacer = PplacerRunner(threads=options.pplacer_threads)  # fix at one thread to keep memory requirements reasonable
        pplacer.run(binFiles, resultsParser, options.output_dir, options.bReducedTree)

        self.timeKeeper.printTimeStamp()

    def treeQA(self, options):
        """QA command"""
        self.logger.info('[CheckM - tree_qa] Assessing phylogenetic markers found in each bin.')

        checkDirExists(options.tree_dir)

        # set HMM file for each bin
        markerSetParser = MarkerSetParser()
        hmmModelInfoFile = os.path.join(options.tree_dir, 'storage', DefaultValues.PHYLO_HMM_MODEL_INFO)
        binIdToModels = markerSetParser.loadBinModels(hmmModelInfoFile)

        # calculate marker gene statistics
        RP = ResultsParser(binIdToModels)
        binStats = RP.analyseResults(options.tree_dir,
                                          DefaultValues.BIN_STATS_PHYLO_OUT,
                                          DefaultValues.HMMER_TABLE_PHYLO_OUT)

        # determine taxonomy of each bin
        
        treeParser = TreeParser()
        treeParser.printSummary(options.out_format, options.tree_dir, RP, options.bTabTable, options.file, binStats)

        if options.file != '':
            self.logger.info('QA information written to: ' + options.file)

        self.timeKeeper.printTimeStamp()

    def lineageSet(self, options, db=None):
        """Lineage set command"""
        self.logger.info('[CheckM - lineage_set] Inferring lineage-specific marker sets.')

        checkDirExists(options.tree_dir)

        # set HMM file for each bin
        markerSetParser = MarkerSetParser()
        hmmModelInfoFile = os.path.join(options.tree_dir, 'storage', DefaultValues.PHYLO_HMM_MODEL_INFO)
        binIdToModels = markerSetParser.loadBinModels(hmmModelInfoFile)

        # calculate marker gene statistics

        resultsParser = ResultsParser(binIdToModels)
