def cigarTuplesToOperations(tuples, includeHardClip=True):
    """
    Produce a sequence of CIGAR operations given a cigar tuples list
    from a pysam AlignedRead instance.

    See https://pysam.readthedocs.io/en/latest/
        api.html#pysam.AlignedSegment.cigartuples

    @param tuples: A C{list} of 2-tuples, each containing an operation
        and an C{int} count.
    @param includeHardClip: If C{True}, yield hard clipping operations,
        else elide them.
    """
    for operation, count in tuples:
        if operation != CHARD_CLIP or includeHardClip:
            for _ in range(count):
                yield operation def softClippedOffset(offset, pairs, cigarOperations):
    """
    Get the offset in the reference of a base that has been soft-clipped
    in a query.

    @param offset: The offset (in C{pairs}) of the soft-clipped base (which
       has a to-be-determined offset in the reference).
    @param pairs: A C{list} of C{int} (queryOffset, referenceOffset) pairs
       for a pysam C{pysam.AlignedSegment}, as returned by
       C{read.get_aligned_pairs} where the read is returned by the pysam
       C{fetch} method.
    @param cigarOperations: A C{list} of pysam CIGAR operations.
    @return: The C{int} reference offset of the soft-clipped base.
    """
    assert cigarOperations[offset] == CSOFT_CLIP
    assert pairs[offset][1] is None

    # Look back.
    count = 0
    for pair, cigarOperation in zip(pairs[offset::-1],
                                    cigarOperations[offset::-1]):
        if cigarOperation == CSOFT_CLIP:
            count += 1
        else:
            _, referenceOffset = pair
            # If we have a reference offset, then we've found our
            # answer. Otherwise, this has to be a hard-clipped base.
            if referenceOffset is None:
                assert cigarOperation == CHARD_CLIP
                # We could break but won't, as explained above.
            else:
                return referenceOffset + count

    # Look ahead.
    count = 0
    for pair, cigarOperation in zip(pairs[offset:], cigarOperations[offset:]):
        if cigarOperation == CSOFT_CLIP:
            count += 1
        else:
            _, referenceOffset = pair
            # If we have a reference offset, then we've found our
            # answer. Otherwise, this has to be a hard-clipped base.
            if referenceOffset is None:
                assert cigarOperation == CHARD_CLIP
                # We could break here, but if we keep going we can test
                # that the whole rest of the CIGAR is hard clips and will
                # fail if not. That's not needed in a world of perfect code.
            else:
                return referenceOffset - count

    # This should be impossible.
    raise ValueError('Soft-clipped base with no following or preceding '
                     'non-hard-clipped bases.') def insertionOffset(offset, pairs, cigarOperations):
    """
    Get the insertion offset in the reference for a base that is part of an
    insertion.

    @param offset: The offset (in C{pairs}) of the soft-clipped base (which
       has a to-be-determined offset in the reference).
    @param pairs: A C{list} of C{int} (queryOffset, referenceOffset) pairs
       for a pysam C{pysam.AlignedSegment}, as returned by
       C{read.get_aligned_pairs} where the read is returned by the pysam
       C{fetch} method.
    @param cigarOperations: A C{list} of pysam CIGAR operations.
    @return: A 2-tuple containing a C{bool} and the C{int} reference offset of
        the base after the insertion. The C{bool} indicates whether we needed
        to look back in the sequence to find a reference base.
    """
    assert cigarOperations[offset] == CINS
    assert pairs[offset][1] is None

    # Look back.
    for pair, cigarOperation in zip(pairs[offset::-1],
                                    cigarOperations[offset::-1]):
        if cigarOperation != CINS:
            _, referenceOffset = pair
            # If we have a reference offset, then we've found our
            # answer. Otherwise, check as above.
            if referenceOffset is None:
                assert cigarOperation not in CONSUMES_REFERENCE
            else:
                return True, referenceOffset + 1

    # Look ahead.
    for pair, cigarOperation in zip(pairs[offset:], cigarOperations[offset:]):
        if cigarOperation != CINS:
            _, referenceOffset = pair
            # If we have a reference offset, then we've found our
            # answer. Otherwise, check that this is something that does not
            # consume reference bases.
            if referenceOffset is None:
                assert cigarOperation not in CONSUMES_REFERENCE
            else:
                return False, referenceOffset

    # This should be impossible.
    raise ValueError('Inserted base with no following or preceding '
                     'reference bases.') class Bases:
    """
    Manage a collection of (base, quality) pairs for a genome site.
    """
    __slots__ = ('count', 'counts')

    def __init__(self):
        self.count = 0
        self.counts = dict.fromkeys('ACGT', 0)

    def __eq__(self, other):
        return self.count == other.count and self.counts == other.counts

    def __str__(self):
        return f'<Bases count={self.count}, bases={self.counts}'

    __repr__ = __str__

    def __getitem__(self, key):
        return self.counts[key]

    def __len__(self):
        return self.count

    def __add__(self, other):
        new = Bases()
        new.count = self.count + other.count
        for counts in self.counts, other.counts:
            for key, count in counts.items():
                new.counts[key] += count
        return new

    def __iadd__(self, other):
        self.count += other.count
        for key, count in other.counts.items():
            self.counts[key] += count
        return self

    def append(self, base, quality):
        """
        Append a (base, quality) pair.

        @param base: A C{str} nucleotide base.
        @param quality: An C{int} nucleotide quality score.
        @return: C{self}.
        """
        if base != 'N':
            self.count += 1
            self.counts[base] += quality

        return self

    def consensus(self, threshold, minCoverage, lowCoverage, noCoverage):
        """
        Get the base that can be used as part of a consensus.

        If there are sufficient reads, this is the least-ambiguous nucleotide
        code for our bases, given a required homogeneity threshold. Otherwise,
        the low coverage value.

        @param threshold: A C{float} threshold. This fraction, at least, of the
            most-common nucleotides at a site are used to determine the
            consensus nucleotide (or ambiguous symbol if more than one
            nucleotide is required to achieve this threshold). If there is a
            tie in nucleotide counts at a site that causes the threshold to be
            met, all nucleotides of equeal frequncy will be included in the
            ambiguous symbol for that site. This is perhaps better explained
            with an example. See
            https://assets.geneious.com/manual/2020.1/static/GeneiousManualse43.html
            and the corresponding testGeneiousExamplesTie test in
            test/test_dna.py
        @param minCoverage: An C{int} minimum number of reads that must cover a
            site for a consensus base to be called. If fewer reads cover a
            site, the C{lowCoverage} value is used.
        @param lowCoverage: A C{str} indicating what base to use when
            0 < N < minCoverage reads cover a site.
        @param noCoverage: A C{str} indicating what base to use when
            no reads cover the site.
        @return: A C{str} nucleotide code. This will be an ambiguous code if
            the homogeneity C{threshold} is not met.
        """
        return (noCoverage if self.count == 0 else
                (lowCoverage if self.count < minCoverage else
                 leastAmbiguousFromCounts(self.counts, threshold))) @contextmanager
def maybeProgressBar(show, maxValue, prefix):
    """
    A context manager to maybe show a progress bar.

    @param show: If C{True}, yield a progress bar, else a Class with an
        C{update} method that does nothing.
    @param maxValue: The C{int} number of tasks to show progress for.
    @param prefix: A C{str} prefix, to appear at the start of the progress bar.
    """
    if show and os.isatty(2):
        widgets = [
            progressbar.SimpleProgress(format='%(value_s)s/%(max_value_s)s'),
            progressbar.Percentage(format=' %(percentage)3d%%'),
            ' ',
            progressbar.Bar(marker='\x1b[33m#\x1b[39m'),
            ' ',
            progressbar.Timer(format='Elapsed: %(elapsed)s'),
            ' ',
            progressbar.ETA(format='ETA: %(eta)8s'),
        ]
        with progressbar.ProgressBar(max_value=maxValue, widgets=widgets,
                                     prefix=prefix) as bar:
            yield bar
    else:
        class Bar:
            update = staticmethod(lambda _: None)
        yield Bar class DNARead(_NucleotideRead):
    """
    Hold information and methods to work with DNA reads.
    """
    ALPHABET = set('ATCG')

    COMPLEMENT_TABLE = _makeComplementTable(ambiguous_dna_complement) @contextmanager
def samfile(filename):
    """
    A context manager to open and close a SAM/BAM file.

    @param filename: A C{str} file name to open.
    """
    f = AlignmentFile(filename)
    yield f
    f.close() def getReferenceInfo(bam, bamFilename, bamId=None, referenceFasta=None,
                     fastaId=None, quiet=False):
    """
    Get the id and length of a BAM file reference seqeunce.

    @param bam: An open BAM file.
    @param bamFilename: The C{str} file name of the BAM file.
    @param bamId: A C{str} BAM file reference name indicating which aligned
        reads to make a consensus from. If not given, will be inferred
        from the BAM file header.
    @param referenceFasta: A C{str} file name containing the sequence that was
        aligned to in making the BAM file.
    @param fastaId: A C{str} reference name indicating which sequence in
        C{referenceFasta} to use as a reference. Only considered if
        C{referenceFasta} is given. If not given and C{referenceFasta} is,
        the reference id will be inferred from reference names in the BAM
        header, or will be taken as the id of the first sequence in
        C{referenceFasta}.
    @param quiet: If C{True}, suppress diagnostic output.
    @raise UnknownReference: For an unknown reference.
    @raise UnspecifiedReference: If a reference is not given and cannot be
        inferred.
    @raise ReferenceNameMismatchError: If the names of the reference in the BAM
        file and the FASTA reference file (if given) do not match (unless both
        are explicitly given, in which case mismatching names are allowed as
        long as the sequence lengths are the same).
    @raise UnequalReferenceLengthError: If the lengths of the reference in the
        BAM file and the FASTA reference file (if given) do not match
    @return: A 2-tuple of C{str}, C{int} giving the refrence sequence id and
        length.
    """
    bamReferences = set(samReferences(bam))
    checkName = checkLength = True
    inferredBamId = bamId is None
    inferredFastaId = fastaId is None
    referenceLength = None

    if bamId is None:
        inferredBamId = True
    else:
        if bamId not in bamReferences:
            raise UnknownReference(
                f'BAM file {str(bamFilename)!r} does not mention a reference '
                f'with id {bamId!r}. Known references are: '
                f'{", ".join(sorted(bamReferences))}.')

        inferredBamId = False
        tid = bam.get_tid(bamId)
        referenceLength = bam.lengths[tid]

        if fastaId is not None:
            # Both ids have been given explicitly, so don't check that they
            # are identical, assume the user knows what they are doing (we
            # will just check the lengths, if necessary).
            checkName = False

    if referenceFasta is None:
        reference = None
    else:
        if fastaId is None:
            # We're given a reference FASTA file, but no id to use. Look at
            # the sequences in the FASTA and take the first one that has a
            # name and length matching a BAM reference (we could actually
            # look to see if there are more than one and, if so, give an
            # error if their sequences are not identical).
            firstRead = None
            for read in FastaReads(referenceFasta):
                if firstRead is None:
                    firstRead = read
                if read.id in bamReferences:
                    tid = bam.get_tid(read.id)
                    if len(read) == bam.lengths[tid]:
                        reference = read
                        checkLength = checkName = False
                        break
            else:
                # Try the first sequence in the FASTA file.
                reference = firstRead

            if firstRead is None:
                raise UnknownReference(
                    f'The FASTA reference file {str(referenceFasta)!r} '
                    f'contained no sequences.')
        else:
            for read in FastaReads(referenceFasta):
                if read.id == fastaId:
                    reference = read
                    break
            else:
                raise UnknownReference(f'No sequence with id {fastaId!r} '
                                       f'found in {str(referenceFasta)!r}.')

        # Set reference length if it has not already been set due to a
        # given bamId.
        if referenceLength is None:
            referenceLength = len(reference)

    if bamId is None:
        if len(bamReferences) == 1:
            # This is the only possibility.
            bamId = tuple(bamReferences)[0]
            tid = bam.get_tid(bamId)
            referenceLength = bam.lengths[tid]
        elif reference and reference.id in bamReferences:
            tid = bam.get_tid(reference.id)
            if len(reference) == bam.lengths[tid]:
                bamId = reference.id
                checkLength = checkName = False
        else:
            raise UnspecifiedReference(
                f'Could not infer a BAM reference. Available references are: '
                f'{", ".join(sorted(bamReferences))}.')

    if inferredBamId and not quiet:
        print(f'BAM reference id {bamId!r} inferred from context.',
              file=sys.stderr)

    if reference:
        if inferredFastaId and not quiet:
            print(f'FASTA reference id {reference.id!r} inferred from '
                  f'context.', file=sys.stderr)

        if checkName and reference.id != bamId:
            raise ReferenceNameMismatchError(
                f'Reference FASTA sequence name {reference.id!r} does '
                f'not match the BAM id {bamId!r}.')

        if checkLength and len(reference) != referenceLength:
            raise UnequalReferenceLengthError(
                f'Reference FASTA sequence {reference.id!r} has length '
                f'{len(reference)}, but the BAM reference {bamId!r} has '
                f'length {referenceLength}.')

    return bamId, reference, referenceLength class UnspecifiedReference(SamError):
    "Reference sequence not specified." CONSUMES_REFERENCE = {CMATCH, CDEL, CREF_SKIP, CEQUAL, CDIFF} def pct(a, b):
    """
    Format a string showing two integers and what percentage the first
    is of the second.

    @param a: An C{int}, the numerator.
    @param b: An C{int}, the denominator.
    """
    assert 0 <= a <= b
    if b:
        return ('%d/%d (%.3f%%)' %
                (a, b, (a / b if b else 0.0) * 100.0))
    else:
        return '0/0 (0.000%)' @contextmanager
def openOr(filename, mode='r', defaultFp=None, specialCaseHyphen=True):
    """
    A context manager to either open a file or use a pre-opened default file.

    @param filename: If not C{None}, this is the argument to pass to C{open}
        along with C{mode}. If C{None}, C{fp} is used.
    @param mode: The C{str} file opening mode, used if C{filename} is not
        C{None}.
    @param defaultFp: An open file-like object to yield if C{filename} is
        C{None}.
    @param specialCaseHyphen: If C{True}, treat '-' as a C{None} filename
        and yield the C{defaultFp}.
    """
    if filename is None or filename == '-' and specialCaseHyphen:
        yield defaultFp
    else:
        with open(filename, mode) as fp:
            yield fp
import sys
from collections import defaultdict
from pysam import CINS, CDEL, CSOFT_CLIP
from dark.cigar import (
    cigarTuplesToOperations, softClippedOffset, insertionOffset)
from dark.dna import Bases
from dark.progress import maybeProgressBar
from dark.reads import DNARead
from dark.sam import (
    samfile, getReferenceInfo, UnspecifiedReference, CONSUMES_REFERENCE)
from dark.utils import pct, openOr
            if offset in otherBases and (
                    offset not in offsetBases or not offsetBases[offset]):
                result.append(deletionSymbol if otherBases[offset] is None
                              else otherBases[offset])
            else:
                assert offset in offsetBases, (
                    f'Offset {offset} not found in offsetBases or otherBases.')
                bases = offsetBases[offset]
                result.append(bases.consensus(
                    threshold, minCoverage, lowCoverageBase, noCoverageBase))

            bar.update(barCount)

    return ''.join(result)


def getConsensusId(bamId, idLambda):
    """
    Make an id for the consensus sequence.

    @param bamId: A C{str} reference sequence id from the BAM file.
    @param idLambda: A one-argument function taking and returning a sequence
        id. This can be used to set the id of the consensus sequence based
        on the id of the reference sequence. The function will be called with
        the id of the BAM reference sequence.
    @return: A C{str} sequence id for the consensus.
    """
    if idLambda is None:
        consensusId = f'{bamId}-consensus'
    else:
        idLambda = eval(idLambda)
        consensusId = idLambda(bamId)

    return consensusId


def consensusFromBAM(
        bamFilename, bamId=None, referenceFasta=None, fastaId=None,
        consensusId=None, idLambda=None, threshold=0.8, minCoverage=1,
        lowCoverage='n', noCoverage='N', deletionSymbol='-',
        deletionThreshold=0.5, ignoreQuality=False, insertionCountThreshold=5,
        strategy='fetch', includeSoftClipped=False, compareWithPileupFile=None,
        progress=False, quiet=False):
    """
    Build a consensus sequence from a BAM file.

    @param bamFilename: the BAM file.
    @param bamId: A C{str} BAM file reference name indicating which aligned
        reads to make a consensus from. If not given, will be inferred
        from the BAM file header.
    @param referenceFasta: A C{str} file name containing the sequence that was
        aligned to in making the BAM file.
    @param fastaId: A C{str} reference name indicating which sequence in
        C{referenceFasta} to use as a reference. Only considered if
        C{referenceFasta} is given. If not given and C{referenceFasta} is,
        the reference id will be inferred from reference names in the BAM
        header, or will be taken as the id of the first sequence in
        C{referenceFasta}.
    @param consensusId: The C{str} id to use in the consensus sequence. If not
        given, the BAM reference id with '-consensus' appended will be used.
    @param idLambda: A one-argument function taking and returning a sequence
        id. This can be used to set the id of the consensus sequence based
        on the id of the reference sequence. The function will be called with
        the id of the BAM reference sequence.
    @param threshold: A C{float} threshold. This fraction, at least, of the
        most-common nucleotides at a site are used to determine the consensus
        nucleotide (or ambiguous symbol if more than one nucleotide is
        required to achieve this threshold). If there is a tie in nucleotide
        counts at a site that causes the threshold to be met, all nucleotides
        of equeal frequncy will be included in the ambiguous symbol for that
        site. This is perhaps better explained with an example. See
        https://assets.geneious.com/manual/2020.1/static/GeneiousManualse43.html
        and the corresponding testGeneiousExamplesTie test in test/test_dna.py
    @param minCoverage: An C{int} minimum number of reads that must cover a
        site for a consensus base to be called. If zero reads cover a site, the
        C{noCoverage} value is used or if the number is greater than zero but
        less than C{minCoverage}, the C{lowCoverage} value is used.
    @param lowCoverage: A C{str} indicating what to do when some reads cover a
        site, but fewer than C{minCoverage}. Either 'reference' or a single
        character (e.g., 'N').
    @param noCoverage: A C{str} indicating what to do when no reads cover a
        reference base. Either 'reference' or a single character (e.g., 'N').
    @parm deletionSymbol: The C{str} to insert in the consensus when a deleted
        site is detected.
    @param deletionThreshold: If some reads have a deletion at a site and some
        do not, call the site as a deletion if C{float} the fraction of reads
        with the deletion is at least this value.
    @param insertionCountThreshold: The C{int} number of reads that must have
        an insertion at an offset in order for the insertion to be called in
        the consensus.
    @param ignoreQuality: If C{True}, ignore quality scores.
    @param strategy: A C{str} consensus-making strategy.
    @param includeSoftClipped: Include information from read bases that were
        marked as soft-clipped by the algorithm that made the BAM file.
    @param compareWithPileupFile: If C{True}, compare the base counts from the
        pysam fetch method with those of the pileup methods. This pays no
        attention to insertions. A summary of the result is written to this
        file.
    @param progress: If C{True}, display a progress bar on standard error.
    @param quiet: If C{True}, suppress diagnostic output. Note that this will
        silence warnings about differing reference names.
    @raise UnspecifiedReference: If no id is provided to indicate which BAM
        file reference to call a consensus for.
    @raise UnknownReference: If a requested reference id is unknown.
    @raise UnequalReferenceLengthError: If the passed reference does not have a
        length identical to the length mentioned in the BAM file.
    @raise ReferenceNameMismatchError: If the name of the FASTA reference
        sequence and the BAM reference do not agree (this is not raised if both
        ids are given explicitly).
    @return: A C{Read} instance with the consensus sequence.
    """
    if referenceFasta is None:
        if lowCoverage == 'reference':
            raise UnspecifiedReference('lowCoverage is "reference" but no '
                                       'reference FASTA file was given.')

        if noCoverage == 'reference':
            raise UnspecifiedReference('noCoverage is "reference" but no '
                                       'reference FASTA file was given.')


    with samfile(bamFilename) as bam: